title: PyMC3 performance issues with No U-Turn Sampler (NUTS): less than 2 iterations per second with simple model
body: <p>I am trying to understand what is the problem with the following code:</p>
  
  <pre><code>import pymc3 as pm
  import theano as t
  
  X = t.shared(train_new)
  features = list(map(str, range(train_new.shape[1])))
  with pm.Model() as logistic_model:
      glm = pm.glm.GLM(X, targets, labels=features,
                       intercept=False, family='binomial')
      trace = pm.sample(3000, tune=3000, jobs=-1)
  </code></pre>
  
  <p>The dataset is by no means big: its shape is <code>(891, 13)</code>. Here is what I concluded on my own:</p>
  
  <ul>
  <li>the problem is surely not the hardware because the performance is the same both on my laptop and on a c4.2xlarge AWS instance;  </li>
  <li>it cannot be <code>theano.shared</code> because if I remove it the performance is again the same;</li>
  <li>the problem does not appear to be in <code>pymc3.glm.GLM</code> because when I manually build the model (which is probably simpler than the one in <code>GLM</code>) the performance is just as terrible:
  <br/><br/><pre><code>with pm.Model() as logistic_model:
      invlogit = lambda x: 1 / (1 + pm.math.exp(-x))
      σ = pm.HalfCauchy('σ', beta=2)
      β = pm.Normal('β', 0, sd=σ, shape=X.get_value().shape[1])
      π = invlogit(tt.dot(X, β))
      likelihood = pm.Bernoulli('likelihood', π, observed=targets)
  </pre></code></li>
  </ul>
  
  <p>It starts at around <code>200 it/s</code> and the quickly falls to <code>5 it/s</code>. After half sampling, it decreases further to around <code>2 it/s</code>. This is a serious problem, as the model barely converges with a couple of thousands of samplings. I need to perform many more samples than what this situation currently allows.</p>
  
  <p>This is the log:</p>
  
  <pre><code>Auto-assigning NUTS sampler...
  Initializing NUTS using jitter+adapt_diag...
   99%|█████████▊| 5923/6000 [50:00&lt;00:39,  1.97it/s]
  </code></pre>
  
  <p>I tried with <code>pm.Metropolis()</code> as step, and it was a bit faster but it didn't converge.</p>
  
  <p><strong>MWE</strong>: a file with a minimal working example showing the issue and the data is here:
  <a href="https://gist.github.com/rubik/74ddad91317b4d366d3879e031e03396" rel="nofollow noreferrer">https://gist.github.com/rubik/74ddad91317b4d366d3879e031e03396</a></p>