major slowdown in distribution.random()

[hvasbath]

I abserved a major slowdown in using the distributions.random() method to create an initial starting point population for metropolis chains the following way:

# create initial population
vars = inputvars(vars)
self.population = []
for i in range(self.n_chains):
    dummy = pm.Point({v.name: v.random() for v in vars},
                                                            model=model)
    self.population.append(dummy)

On the most recent version it takes: 260.339263s

And on the commit:
commit 1c9adc6
Author: Thomas Wiecki
Date: Tue Nov 29 11:30:52 2016 +0100

MAINT Re-enable trust_input.

it takes 6.925914 s.

I dont know what you changed where, unfortunately I couldnt keep up with all the commits. Have there been any major changes in the random or subsequently related methods?

The particular distribution that is being called there is the Uniform distribution.

[twiecki]

I bet that's #1585. Can you test that commit? CC @ferrine

[ferrine]

<method 'poll' of 'select.poll' objects> Don't know what does it mean but it takes 74% of runtime in test_distributions_random

[ferrine]

@hvasbath could you please prepare a reproducible example so that I can profile it on my local machine? From the code above it is not obvious what random variables are sampled.

[hvasbath]

I wrote it in the last line in the post above- it is the uniform distribution. Maybe today evening I can extent it a little. But there is really not much more to it- it is n_chains = 1000.

[ferrine]

I cannot reproduce slowdown (branch implement_opvi, for this case is the same as master)

with pm.Model() as model:
    for name in pm.distributions.__all__:
        obj = getattr(pm.distributions, name)
        if isinstance(obj, type) and issubclass(obj, pm.distributions.Distribution):
            try:
                v = obj(name, shape=(10, 10))
                v.random()
                print('testing %s' % name)
                %timeit pm.Point({v.name: v.random()}, model=model)
            except:
                pass

testing Uniform
100 loops, best of 3: 9.92 ms per loop
testing Normal
100 loops, best of 3: 10.4 ms per loop
testing Lognormal
100 loops, best of 3: 4.19 ms per loop
testing HalfNormal
100 loops, best of 3: 9.81 ms per loop
testing SkewNormal
100 loops, best of 3: 16.4 ms per loop
testing Triangular
100 loops, best of 3: 5.82 ms per loop

[hvasbath]

Theano version: '0.9.0rc1.dev-183b5f1a467012b8b15b72b2b22babbb576ee527'

It scales quite linearly with number of variables and apparently also what is inside the model formulation. Now on this simple toy example:

import pymc3 as pm
from pymc3.model import modelcontext
from pymc3.theanof import inputvars
from time import time


n_chains = 1000

with pm.Model() as model:
    var1 = pm.Uniform('v1', shape=(1,), lower=0, upper=10, transform=None)
    var2 = pm.Uniform('v2', shape=(1,), lower=-10, upper=20, transform=None)
    var3 = pm.Uniform('v3', shape=(1,), lower=-20, upper=10, transform=None)
    var4 = pm.Uniform('v4', shape=(1,), lower=0, upper=20, transform=None)
    var5 = pm.Uniform('v5', shape=(1,), lower=0, upper=10, transform=None)
    var6 = pm.Uniform('v6', shape=(1,), lower=-10, upper=200, transform=None)
    var7 = pm.Uniform('v7', shape=(1,), lower=0, upper=10, transform=None)
    var8 = pm.Uniform('v8', shape=(1,), lower=0, upper=20, transform=None)
    llk = var1 + var2 + var3 + var4 + var5 + var6 + var7 + var8

def get_init_pop(model, n_chains):
    model = modelcontext(model)
    vars = inputvars(model.vars)
    print vars
    population = []

    t0 = time()
    for i in range(n_chains):
        dummy = pm.Point({v.name: v.random() for v in vars}, model=model)
        population.append(dummy)

    t1 = time()
    print('time for initialising population', t1 - t0)
    return population

population = get_init_pop(model, n_chains)

On the current master:
[v8, v7, v6, v5, v4, v3, v2, v1]
('time for initialising population', 20.776487827301025)

On the branch mentioned above from November:
[v8, v7, v6, v5, v4, v3, v2, v1]
('time for initialising population', 0.8033387660980225)

[hvasbath]

Any progress on this @ferrine ?

[ferrine]

Sorry. I had a lot of work recently and there was no time to dive deeper. I can check it till weekends

[ferrine]

@hvasbath thank you for providing standalone example. I've done git bisect in master branch. First bad commit is 45c11f0.

-        self.lower = lower
-        self.upper = upper
+        self.lower = lower = tt.as_tensor_variable(lower)
+        self.upper = upper = tt.as_tensor_variable(upper)

After these changes this part of code is executed differently (param now has name attribute) and it can cause slowness. I'm not sure if we can revert the commit. Your thoughts @twiecki?

[twiecki]

@ferrine Good catch. These changes were important to keep compatibility with the theano floatX setting.

Not sure what's causing the slow-down, are we wrapping the parameters in a theano.function call every iteration? Can we get around doing that?

[hvasbath]

@ferrine thanks for taking the time and further investigating!

why do you have to use as_tensor_variable? isnt simply casting to floatX doing the job? Yes there is a theano function call every iteration, but it was like that also in the old version ...

[twiecki]

We could, I think, but I thought it would be cleaner to just have everything be converted to theano objects as early as possible. I still believe that, so we should figure out if we can't remove the bottleneck with theano constants being present.

[hvasbath]

Any progress on this?

[ferrine]

Unfortunately no progress
I have no idea what to do with 45c11f0

[hvasbath]

I am still convinced as stated above that rather casting the values to the respective theano.config.floatX type is enough they dont need to be theano tensors right away. Then the random method is much faster if it works on non-symbolic variables.
What is the reason again for transforming them to theano tensors @twiecki ? I mean they will be still in the wrong float format as as_tensor_variable doesnt affect the dtype.