# Performance issue with Turing (but not Optimization)

**URL:** <https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [June 3, 2025, 8:16pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611 "2025-06-03T20:16:13Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bill\_Holmes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bill_holmes/32/47921_2.png) [@Bill\_Holmes](https://discourse.julialang.org/u/Bill_Holmes)\
**Post date:** [June 3, 2025, 8:16pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/1 "2025-06-03T20:16:13Z")

</div>

Hello all Turing experts. I am having a vexing issue with using NUTS sampling to sample from the posterior of a semi-complicated model involving ODE’s. This is a performance issue that I can’t track down. Unfortunately I can’t describe the whole modeling problem or boil it down to a minimum working example, so here are a few details.

1. I have a biological modeling problem involving a few dozen differential equations for which I am interested in only a few parameters (see below for turing model outline).

2. All of the complexity of the model resides in the custom **logp\_fun** that computes the sum log likelihood of data given parameters.

3. I have been careful in the construction of this function. I have ensured that types can be inferred and that the _remake_ process for continually re-simulating ODE’s follows SciML guidance for optimization contexts. Due to the structure of my data, I am multi-threading within this function (so no MCMCThreads() later). I have however ensured there are no data-races or any threading issues (simple use of Threads.@threads on a loop).

4. I have tested this core function alone on behavior and performance. Computing the function alone is ~1-2ms and computing the gradient is also ~1-2ms (from BenchmarkTools).

5. Optimizing this model using Optimization.jl behaves and performs as expected. Using _Adam_ I get ~3ms / step. Using BBO I get ~1ms / step (no gradient here). And the optimization results themselves are what I expect.

6. To help Turing out, I first use Optimization.jl to get a MLE estimate. I use that as an initial estimate to get a MAP estimate using _maximum\_a\_posteriori_. I then use that MAP estimate as the _initial\_params_ for sampling. So I know that the sampler is starting in the support of the posterior (I also have print checks in _logp\_fun_ to tell me if the model has wandered and it does not).

7. However using NUTS within turing takes 8 sec / sample (single chain, 150 samples in 1293 sec).

To reiterate, the **logp\_fun** and its gradient take 1-3ms to compute. Adam takes about 3ms / step. But NUTS is taking 8sec per sample (1000x longer). I realize that NUTS requires multiple evaluations per sample, but 1000x? Any ideas?

I realize there is not much to work with here since the logp\_fun is too complicated to include here. But the function and its gradients clearly function well both on their own and within a gradient based optimizer (Adam). So there must be something about the NUTS algorithm or the Turing implementation of it that is leading to extremely slow sampling.

```julia
@model function construct_model(logp_fun)

    # Sampling logs to work on real line. Parameter back-conversion in the logp_fun.
    log_p1 ~ Normal(-4.5,0.5)
    log_p2 ~ Normal(0.0,1.0)
    log_p3 ~ Normal(0.0,1.0)
    log_p4 ~ Normal(-4.5,0.5)
    log_σ ~ Normal(-0.7,1.0)

    σ = exp(log_σ)

    params_vec = [log_p1, log_p2, log_p3, log_p4, σ];

     # logp_fun is hand crafted and cannot be reduced to a minimum working example.
    Turing.@addlogprob! logp_fun(params_vec)
end

model = construct_model(logp_fun); #logp_fun defined elsewhere

# Optimization.jl used to get initial MLE estimate. Global optim (BBO) required so Turing not used for this.
map_est = maximum_a_posteriori(model, 
    initial_params = MLE_optim.u,
);

chain_NUTS = sample(model, NUTS(0.7), 100; initial_params=map_est.values.array)

```

---

<div class="post-metadata">

**Author:** ![rto](https://avatars.discourse-cdn.com/v4/letter/r/e480ec/32.png) [@rto](https://discourse.julialang.org/u/rto)\
**Post date:** [June 4, 2025, 5:51am UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/2 "2025-06-04T05:51:14Z")

</div>

Some ideas:

- Initialising at the MAP can often make things worse ([https://statmodeling.stat.columbia.edu/2024/05/24/hmc-fails-when-you-initialize-at-the-mode/](https://statmodeling.stat.columbia.edu/2024/05/24/hmc-fails-when-you-initialize-at-the-mode/)). It is worth trying without specifying any initial params
- You could try sampling with HMC with an arbitrary step size so that you aren’t including the NUTS warmup steps in your timing
- @profview will show where things are slowing down internally

---

<div class="post-metadata">

**Author:** ![scheidan1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/scheidan1/32/24771_2.png) [@scheidan1](https://discourse.julialang.org/u/scheidan1)\
**Post date:** [June 4, 2025, 8:32am UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/3 "2025-06-04T08:32:15Z")

</div>

We had some troubles with the combination of NUTS and ODE solvers that use adaptive time stepping. The adaptive time-stepping can result in “noisy” gradients which then forces NUTS to make tiny steps.

I do not know is there is a formal investigation of this unfortunate interactions. A bit related is this paper, although they focus on optimization and do not use automatic differentiation:

Kavetski, D., Clark, M.P., 2010. Ancient numerical daemons of conceptual hydrological modeling: 2. Impact of time stepping schemes on model analysis and prediction. Water Resources Research 46. [https://doi.org/10.1029/2009WR008896](https://doi.org/10.1029/2009WR008896)

---

<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:** [June 4, 2025, 8:35am UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/4 "2025-06-04T08:35:43Z")

</div>

The noise is related to the tolerance of the ODE solver. So generally for inverse problems the tolerance should be set to something like `abstol=1e-8,reltol=1e-8` or lower. Did you do that in your case?

(Generally if this is your issue though, you would see the same issue with BFGS/LBFGS since those also do not work well with stochasticity in the gradient, but from OP it seems that might not apply here)

---

<div class="post-metadata">

**Author:** ![scheidan1](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/scheidan1/32/24771_2.png) [@scheidan1](https://discourse.julialang.org/u/scheidan1)\
**Post date:** [June 4, 2025, 9:22am UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/5 "2025-06-04T09:22:39Z")

</div>

The noise definitively depends on the tolerances. We did some experiments and I remember that the noise of the gradient was also depending on the solver algorithm (I would need to look up the details, it a while ago)

At least for our toy model, deactivating the adaptive time stepping has sometime lead to bias, but we got always a smooth gradient even for low tolerances.

---

<div class="post-metadata">

**Author:** ![Bill\_Holmes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bill_holmes/32/47921_2.png) [@Bill\_Holmes](https://discourse.julialang.org/u/Bill_Holmes)\
**Post date:** [June 4, 2025, 4:47pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/6 "2025-06-04T16:47:10Z")

</div>

Thanks for the comments and suggestions. I’m hearing two potential issues.

1. Initializing at the MAP estimate may lead to issues either I) due to starting at highest density point or II) poor adaptation phase since gradients should be small at the MAP. I don’t think this is the case. I attempted to move the initial\_params a bit from the MAP and things did not improve. The case @rto mentioned is for high dimensional problems. This is 5-dim.

2. @ChrisRackauckas suggested lower ODE tolerances. This appears to work, but with a strange quirk. I tried three cases (all using AutoTsit5(Rosenbrock23())).  
i) abstol = reltol = 1e-8. The logpfun (and gradient) take ~10ms. NUTS takes ~ **180ms / sample**.  
ii) abstol = reltol = 1e-4. The logpfun and its gradient (ForwardDiff) take ~1.5ms. NUTS takes ~ **180ms / sample**.  
iii) No tolerences included. In this case, the target function and its gradient take ~1.5ms. NUTS takes ~ **8sec / sample**. Further, NUTS misbehaves in the sense that its progress seems to change erratically (based on the progress printing)

So here are the takehomes (for this specific problem)

2a) Specifying tolerances indeed improves NUTS sampling. However, tight tolerances aren’t needed for sampling to proceed. The act of specifying them seems to have some effect.

2b) Tight tolerances indeed slow the logdensity and its gradient evaluation (factor of 7). And the use of Optimization.jl for a MLE estimate slows accordingly (not shown). However, NUTS sampling takes about the same time.

2c) Posteriors obtained with tighter and looser tolerances are near identical, though the ESS is a bit better for the tighter tolerance. The posterior (has to run overnight!!!) without specifying tolerances is terrible.

Not sure the source of this. And the testing was not exhaustive. But it appears specifying tolerances in the ODE solve _may be_ important when using NUTS, though the results aren’t that sensitive to the tolerances.

Would be great if anyone has insights into this. But this is good enough for now.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [June 4, 2025, 5:01pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/7 "2025-06-04T17:01:40Z")

</div>

> [@Bill\_Holmes](#):
>
> To reiterate, the **logp\_fun** and its gradient take 1-3ms to compute. Adam takes about 3ms / step. But NUTS is taking 8sec per sample (1000x longer). I realize that NUTS requires multiple evaluations per sample, but 1000x? Any ideas?

For a complicated posterior NUTS can need quite a lot of evaluations per sample. Turing stores some internal information from the sampler in the chain. You might want to check in particular, `chain_NUTS[info]` where `info ∈ [:n_steps, :acceptance_rate, :tree_depth, :step_size]` to get some additional insights and diagnose problematic sampling.

---

<div class="post-metadata">

**Author:** ![Bill\_Holmes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bill_holmes/32/47921_2.png) [@Bill\_Holmes](https://discourse.julialang.org/u/Bill_Holmes)\
**Post date:** [June 4, 2025, 7:05pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/8 "2025-06-04T19:05:15Z")

</div>

Thanks for the suggestion @bertschi. I wasn’t sure how to access this. I looked at these metrics for the case where the ODE is solved 1) with and 2) without the included tolerance.

The first image is with **abs\_tol = rel\_tol = 1e-4**. The second is with **no tolerance** provided. As expected by NUTS grinding to a halt, each sample requires an enormous number of steps and a deep tree when no tolerance is included. The _n\_steps_, _tree\_depth_ and _step\_size_ are comparable for a tolerance of 1e-8.

The question is why? A tolerance of 1e-4 is not exactly stringent.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/4/7/4741af99175b5fa1148b8df5853a100f211fa9df.png)

 ![image](https://global.discourse-cdn.com/julialang/original/3X/3/6/3698fe5907efd7774067cae26343ca60b527945f.png)

---

<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:** [June 4, 2025, 7:17pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/9 "2025-06-04T19:17:39Z")

</div>

It’s more stringent than the defualt reltol=1e-3. I say for inverse problems a safe default is 1e-8 since it just needs to be “low enough”: the defaults are low enough for forward problems and a basic plot but not for inverse problems.

---

<div class="post-metadata">

**Author:** ![Bill\_Holmes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bill_holmes/32/47921_2.png) [@Bill\_Holmes](https://discourse.julialang.org/u/Bill_Holmes)\
**Post date:** [June 4, 2025, 7:35pm UTC](https://discourse.julialang.org/t/performance-issue-with-turing-but-not-optimization/129611/10 "2025-06-04T19:35:14Z")

</div>

I just tried things with explicit default tolerances of 1e-3 and got the same behavior as with no specified tolerance. Thanks.

For my problem, 1e-4 is sufficient, but NUTS time / sample is the same for 1e-4 and 1e-8 (despite direct simulation time being 7x longer in the 1e-8 case) with slightly better ESS with tighter tolerance. Will of course be problem specific.