Проблемы с производительностью PyMC3 при использовании N-образного сэмплера: менее 2 итераций в секунду с простой моделью

Question

Я пытаюсь понять, в чем проблема со следующим кодом:

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)

Набор данных ни в коем случае не большой: его форма (891, 13), Вот что я сделал самостоятельно:

• проблема, конечно, не в аппаратном обеспечении, потому что производительность одинакова как на моем ноутбуке, так и на экземпляре AWS c4.2xlarge;
• не может быть theano.shared потому что, если я уберу его, производительность снова будет такой же;
• проблема, похоже, не в pymc3.glm.GLM потому что, когда я вручную строю модель (которая, вероятно, проще, чем в GLM) производительность такая же ужасная
  
  

  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)
  

Это начинается около 200 it/s и быстро падает 5 it/s, После половины выборки оно уменьшается до примерно 2 it/s, Это серьезная проблема, так как модель едва сходится с парой тысяч выборок. Мне нужно выполнить гораздо больше сэмплов, чем позволяет эта ситуация.

Это журнал:

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
 99%|█████████▊| 5923/6000 [50:00<00:39,  1.97it/s]

Я пробовал с pm.Metropolis() как шаг, и это было немного быстрее, но не сходилось.

MWE: файл с минимальным рабочим примером, показывающим проблему, и данные здесь: https://gist.github.com/rubik/74ddad91317b4d366d3879e031e03396

Answer 1

Нецентрированная версия модели должна работать намного лучше:

β_raw = pm.Normal('β_raw', 0, sd=1, shape=X.get_value().shape[1])
β = pm.Deterministic('β', β_raw * σ)

Обычно ваш первый импульс, если эффективный размер выборки невелик, должен заключаться не просто в увеличении количества выборок, а в том, чтобы немного поиграть с параметризацией.

Также вы можете использовать tt.nnet.sigmoid вместо вашего собственного invlogit, это может быть быстрее / стабильнее.