Metropolis taking too long to initialize

Questions
1

Hi,

I’m trying to use metropolis method for sampling. It takes too long to initialize the metropolis method.

The code :
step = pm.Metropolis()

This single line takes around 3 minutes to execute.
The sampling itself is pretty fast using metropolis but the initialization is kind of overshadowing all the speed up I gained over slice method.
Is this a known behavior? And can it be improved? In general how long should it takes ideally for the metropolis method to initialize?

I first posted this in Github issues but was suggested to move this here.
Also, shortly I’ll provide the output of code profiling

Thanks

2

Also, more information about the data and the model you are using would be helpful too

3

PFB the profiler output :

********************************************************
 ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001  134.514  134.514 test6.py:1(<module>)
        1    0.000    0.000  118.053  118.053 test6.py:13(run)
       36    0.002    0.000   68.688    1.908 function.py:76(function)
       36    0.003    0.000   68.680    1.908 pfunc.py:283(pfunc)
        1    0.000    0.000   67.645   67.645 arraystep.py:32(__new__)
       32    0.001    0.000   67.645    2.114 metropolis.py:90(__init__)
       32    0.001    0.000   67.477    2.109 metropolis.py:474(delta_logp)
       36    0.003    0.000   64.647    1.796 function_module.py:1722(orig_function)
       36    0.006    0.000   61.128    1.698 function_module.py:1391(__init__)
6524/1261    0.031    0.000   50.851    0.040 opt.py:75(optimize)
       36    0.000    0.000   50.447    1.401 opt.py:92(__call__)
   144/36    0.013    0.000   50.446    1.401 opt.py:213(apply)
32450/32159    0.245    0.000   26.275    0.001 op.py:583(__call__)
     2545    0.004    0.000   26.169    0.010 memoize.py:11(memoizer)
    15893    0.091    0.000   24.787    0.002 op.py:892(make_thunk)
        1    0.000    0.000   24.539   24.539 starting.py:20(find_MAP)
*********************************************************************************

I first create a linear model from formulas, and then run a sampler on it :

glm.GLM.from_formula(formula, data, priors=dict1)
start = find_MAP(fmin=opt.fmin_powell)
step = pm.Metropolis()
trace = sample(n, step, is_cov=True, start=start, njobs=1)

Also, when I try and set njobs>1, it starts using all the cores but the iteration speed decreases significantly from ~22 iterations/sec to ~30 seconds/iteration ans keeps deteriorating with increase in njobs.
Sorry, but passing on data would not be possible I guess :neutral_face:

4

The slowness usually from the find_MAP in my experience. Could you please try the default (NUTS with ADVI init):

glm.GLM.from_formula(formula, data, priors=dict1)
trace = sample(n, njobs=2) #multiple trace for computing Rhat
5

I think I had tried NUTS with ADVI init earlier but the iterations were not as fast as metropolis. But still, I’ll try and run it again and update the results.

6

The faster sample in Metropolis is sometimes a lie :wink: it’s the effective sample that counts. You can have much more effective sample using HMC which means you can use a shorter trace.

7

If I recall, the iterations per second with multiple jobs is only reporting on a single chain (and, you really care about effective samples per second!)