Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to push the for-loop down to numpy

I have the following piece of code doing exactly what i want (it is part of a kriging method). But the problem is that it goes too slow, and i wish to know if there is any option to push the for-loop down to numpy? If i push out the numpy.sum, and use the axis argument there, it speeds up a little bit, but apparently that is not the bottleneck. Any ideas on how i can push down the forloop to numpy to speed it up, or other ways to speed it up?)

# n = 2116
print GRZVV.shape  # (16309, 2116)
print GinvVV.shape  # (2117, 2117) 
VVg = numpy.empty((GRZVV.shape[0]))

for k in xrange(GRZVV.shape[0]):
    GRVV = numpy.empty((n+1, 1))
    GRVV[n, 0] = 1
    GRVV[:n, 0] = GRZVV[k, :]
    EVV = numpy.array(GinvVV * GRVV)  # GinvVV is numpy.matrix
    VVg[k] = numpy.sum(EVV[:n, 0] * VV)

I posted the dimensions of the ndarrays n matrix to clear some stuff out

edit: shape of VV is 2116

like image 534
usethedeathstar Avatar asked Sep 23 '13 11:09

usethedeathstar


2 Answers

You could do the following in place of your loop over k (runtime ~3s):

tmp = np.concatenate((GRZVV, np.ones((16309,1),dtype=np.double)), axis=1)
EVV1 = np.dot(GinvVV, tmp.T)
#Changed line below based on *askewchan's* recommendation
VVg1 = np.sum(np.multiply(EVV1[:n,:],VV[:,np.newaxis]), axis=0)
like image 90
Joel Vroom Avatar answered Oct 01 '22 14:10

Joel Vroom


You are basically taking each row of GRZVV, appending a 1 at the end, multiplying it with GinvVV, then adding up all the elements in the vector. If you weren't doing the "append 1" thing, you could do it all with no loops as:

VVg = np.sum(np.dot(GinvVV[:, :-1], GRZVV.T), axis=-1) * VV

or even:

VVg = np.einsum('ij,kj->k', GinvVV[:, :-1], GRZVV) * VV

How do we handle that extra 1? Well, the resulting vector coming from the matrix multiplication would be incremented by the corresponding value in GinvVV[:, -1], and when you add them all, the value will be incremented by np.sum(GinvVV[:, -1]). So we can simply calculate this once and add it to all the items in the return vector:

VVg = (np.einsum('ij,kj->k', GinvVV[:-1, :-1], GRZVV) + np.sum(GinvVV[:-1, -1])) * VV

The above code works if VV is a scalar. If it is an array of shape (n,), then the following will work:

GinvVV = np.asarray(GinvVV)
VVgbis = (np.einsum('ij,kj->k', GinvVV[:-1, :-1]*VV[:, None], GRZVV) +
          np.dot(GinvVV[:-1, -1], VV))
like image 26
Jaime Avatar answered Oct 01 '22 14:10

Jaime