Проблемы производительности PyMC3 с пробоотборником без разворота (NUTS): менее 2 итераций в секунду с простой моделью

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

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


person rubik    schedule 27.10.2017    source источник


Ответы (1)


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

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

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

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

person aseyboldt    schedule 29.10.2017
comment
Благодарю вас! Я где-то читал, что в таких случаях лучше компенсировать параметризацию, но не смог заставить это работать. В вашей версии я вижу мгновенные 70-кратные улучшения с NUTS. - person rubik; 30.10.2017
comment
Привет, у меня проблемы с производительностью при аналогичной настройке большого набора данных, и я не уверен, связано ли это с этим. Зачем вам бета*сигма? - person swmfg; 30.10.2018