How to write this numpy code in vectorized form?

Question

I have the following function in python which I am not able to figure out how to express in vectorized form. For me, cov is a numpy array of shape (2,2), mu is the mean vector with shape (2,) and xtp is of shape (~50000,2). I know scipy provides scipy.stats.multivariate_normal.pdf but I am trying to learn how to write efficient vectorized code. Please

def mvnpdf(xtp, mu, cov):
    temp = np.zeros(xtp.shape[0])
    i = 0
    length = xtp.shape[0]
    const = 1 / ( ((2* np.pi)**(len(mu)/2)) * (np.linalg.det(cov)**(1/2)) )
    inv = np.linalg.inv(cov)
    while i < length:
        x = xtp[i]-mu
        exponent = (-1/2) * (x.dot(inv).dot(x))
        temp[i] =  (const * np.exp(exponent))
        i+=1
    return temp

hpaulj · Accepted Answer · 2015-09-12 00:10:33Z

4

The only tricky part to vectorize is that double .dot. Let's isolate that:

x = xtp - mu  # move this out of the loop
ddot = [i.dot(inv).dot(i) for i in x]
temp = const * np.exp(-0.5 * ddot)

Put that in your code and see if it produces that same thing.

There are several ways of 'vectorizing' a dot. The one I like to try first is einsum. In my tests this is equivalent:

ddot = np.einsum('ij,jk,ik->i',x,inv,x)

I'd suggest trying it to see if it works and speeds things up. And do play around with these calculations with smaller arrays (not the ~50000) in an interactive shell.

I'm testing things with

In [225]: x
Out[225]: 
array([[  0.,   2.],
       [  1.,   3.],
       ...
       [  7.,   9.],
       [  8.,  10.],
       [  9.,  11.]])
In [226]: inv
Out[226]: 
array([[ 1.,  0.],
       [ 0.,  1.]])

Since this is a learning exercise I'll leave the details to you.

With (2,2), the calculations one cov might be faster it done explicitely rather than with the det and inv functions. But it's that length iteration that's the time consumer.

answered Sep 12, 2015 at 0:10

hpaulj

233k14 gold badges260 silver badges392 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

enitihas Over a year ago

Thanks a lot. einsum looks awesome. The resulting code is a lot faster .Any good resources to learn on using enisum?

Alex Riley Over a year ago

@enitihas: there are relatively few resources scattered around the web, but here's a SO question with some answers and links: stackoverflow.com/questions/26089893/understanding-einsum-numpy

Collectives™ on Stack Overflow

How to write this numpy code in vectorized form?

1 Answer 1

2 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related