How to push the for-loop down to numpy

Question

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

Answer 1

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)

Answer 2

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))