Speed up double loop in Python

Question

Is there a way to speed up a double loop that updates its values from the previous iteration?

In code:

def calc(N, m):
    x = 1.0
    y = 2.0
    container = np.zeros((N, 2))
    for i in range(N):
      for j in range(m):
        x=np.random.gamma(3,1.0/(y*y+4))
        y=np.random.normal(1.0/(x+1),1.0/sqrt(x+1))
      container[i, 0] = x
      container[i, 1] = y
    return container

calc(10, 5)

As you can see, the inner loop is updating variables x and y while the outer loop starts with a different value of x each time. I don't think this is vectorizable but maybe there are other possible improvements.

Thanks!

Answer 1

I don't think it's going to add up to any important speed up, but you can save some function calls if you generate all your gamma and normally distributed random values at once.

Gamma functions have a scaling property, so that if you draw a value x from a gamma(k, 1) distribution, then c*x will be a value drawn from a gamma(k, c) distribution. Similarly, with the normal distribution, you can take a y value drawn from a normal(0, 1) distribution and convert it into a value drawn from a normal(m, s) distribution doing x*s + m. So you can rewrite your function as follows:

def calc(N, m):
    x = 1.0
    y = 2.0
    container = np.zeros((N, 2))
    nm = N*m
    gamma_vals = np.random.gamma(3, 1, size=(nm,))
    norm_vals = np.random.normal(0, 1, size=(nm,))
    for i in xrange(N):
        for j in xrange(m):
            ij = i*j
            x = gamma_vals[ij] / (y*y+4)
            y = norm_vals[ij]/np.sqrt(x+1) + 1/(x+1)
        container[i, 0] = x
        container[i, 1] = y
    return container

If the actual parameters of your distributions had a simpler expression, you may actually be able to use some elaborate form of np.cumprod or the like, and spare yourself the loops. I am not able to figure out a way of doing so...

Answer 2

Does this work?

for i in xrange(N): 
   # xrange is an iterator, range makes a new list.
   # You save linear space and `malloc`ing time by doing this

    x += m*y # a simple algebra hack. Compute this line of the loop just once instead of `m` times
    y -= m*x
    y *= -1 # another simple algebra hack. Compute this line of the loop just once instead of `m` times
    container[i,0] = x
    container[i,1] = y
return container