Proposal generation

[hvasbath]

Is horribly slow! That was a main surprise for me, especially for the Multivariate Proposal for a Covariance matrix of 450x450 it took 0,3s!!! For each step!
Fortunately, there is a solution to that. Just the initialisation of the RNG is so slow in Numpy!
You can call the function with a size argument and generate your proposal steps before the loop! That resulted for me in a sampling that is unbelievable 2 times faster! Maybe with the other distributions it is not as significant, but we might want to consider doing that for the other step_methods at least Metropolis as well. I dont know much about the other methods you guys implemented so far.
Depending on your opinions, of course!
An example of how I did it in the ATMCMC, is in these commits:
98ffaeb
b4ff0e1

[twiecki]

Wow, that's significant. I wonder if this were faster if we cached a scipy.distributions object, e.g.:

class MultivariateNormalProposal(Proposal):
    def __init__(self, s):
        self.random_gen = scipy.stats.multivariate_normal(mean=np.zeros(self.s.shape[0]), cov=self.s)

    def __call__(self):
        return self.random_gen.rvs()

In my local experiments, that was way way faster than calling the numpy.random version.

[twiecki]

I get 8.36ms for numpy, and 187µs with caching, so a 44x speed-up.

[hvasbath]

Yes probably thats ok like that. I tested to use the same numpy function but using the size argument and create 200000 samples from that distribution. That took15.9s for 200.000 samples. Resulting in 7.95micro seconds per sample (speedup of 37735x :) compared to 300ms per sample). Of course that becomes more efficient or less efficient with the size increasing/ decreasing. But your version still uses less memory. How about combining the scipy.distribution with the array pre-generation. Shouldnt this even be faster. I will try it.

[hvasbath]

No it doesnt make a difference.
%time out = random_gen.rvs(1000)
CPU times: user 304 ms, sys: 4 ms, total: 308 ms
%time out = np.random.multivariate_normal(mean=np.zeros(450),cov=np.eye(450),size=1000)
CPU times: user 300 ms, sys: 0 ns, total: 300 ms
But again, the key is creating 'draws' samples at once, rather than repeated calling!

[twiecki]

The scaling get's adjusted in e.g. Metropolis so pre-caching won't work there I think. So we can only use it for ATMCMC?

I realized I only used a 1-D distribution in my measurements above. If I use 450 like you, I also don't get a speed-up with precreating the scipy distribution. I wonder if theano random number generation would be any faster in this context...

[springcoil]

That's an interesting question. I've not used Theano random number
generation does anyone like @fonnesbeck have any experience with it?

On Wed, Mar 23, 2016 at 2:15 PM, Thomas Wiecki notifications@github.com
wrote:

The scaling get's adjusted in e.g. Metropolis so pre-caching won't work
there I think. So we can only use it for ATMCMC?

I realized I only used a 1-D distribution in my measurements above. If I
use 450 like you, I also don't get a speed-up with precreating the scipy
distribution. I wonder if theano random number generation would be any
faster in this context...

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
#1034 (comment)

Peadar Coyle
Skype: springcoilarch
www.twitter.com/springcoil
peadarcoyle.wordpress.com

[hvasbath]

I also do the scaling in ATMCMC that is no problem, because you can still multiply the scaling on your draw. See lines 121 and 122 in here: 98ffaeb. Theano random numers are even slower than numpy and maximum as fast as numpy as it uses numpys RNGs ;)
Theano/Theano#1237
That is something they want to improve for the next big 0.9 release.

[twiecki]

OK, nevermind theano then.

I suppose the question is how much this is slowing down metropolis. The work-around you have seems to work quite well for ATMCMC so we could just have it there at first.

[hvasbath]

The point is that it speeds up Metropolis rather than slowing it down.

[twiecki]

oh so you did measure that. The problem with the manual rescaling is that only works for a normal, right?

[twiecki]

not e.g. poisson

[hvasbath]

Sorry I am lost here. Manual rescaling? You mean the adaptive tuning? That works on any proposal.
Any numpy distribution has a size argument which you can use similarly as I did with the MvNormal.
The tuning is not changed/ influenced at all.

[twiecki]

You're right, we do the same thing here (this is what I meant with manual scaling): https://github.com/pymc-devs/pymc3/blob/master/pymc3/step_methods/metropolis.py#L112

I got confused because I thought we pass in the scaling to the proposal but we're doing what you're doing. I thought maybe multiplying might not always be the correct thing to do, e.g. the poisson uses a different scaling effect: https://github.com/pymc-devs/pymc3/blob/master/pymc3/step_methods/metropolis.py#L44 but I think this should work out fine.

Do you want to do a PR with the precaching of proposal draws?

[hvasbath]

Yes I will do that ;) . Practicing my github skills ...