My apologies, I realised that I should’ve used: at.tril(at.ones(timeD.shape), k=-1) , with k=-1. In this case the model gets stuck at the same place as with the numpy method at slightly over 30-40% . I’ve tried running with 1 core only, but the problem persists, get’s stuck at around 30%. I’ve changed the init method to “adapt_diag” but the same problem occurs. And using init=“advi” the sampling gets stuck at around 60% . I’m clueless 
. May be something internal to pytensor or PyMC? I’m not sure whether JIT would help much with this, or maybe it could?