1

I'm looking for a more efficient way to fill a 2d numpy array than a double for-loop. The issue I am having is that the array values are dependent on several other arrays.

In the following code k and d are integers, and y and result are arrays. Any help would be greatly appreciated.

for i in xrange(1,d):
    for j in xrange(k-i-1, k+1):
        result[j][i] = ((x - y[j])/(y[j+i] - y[j]))*result[j,i-1] + ((y[j+i+1]-x)/(y[j+i+1] - y[j+1]))*result[j+1,i-1]
2
  • are you modifying result in place? That makes we worry about results[j+1, i-1]. Commented Oct 20, 2013 at 14:58
  • I am, but result is initially declared as result = numpy.zeros((len(y),d), so I am not accessing anything that has not already been declared. Commented Oct 20, 2013 at 15:02

2 Answers 2

4

You seem to be updating your result array one column at a time, using data from the previously updated column. That makes it hard to vectorize the outer loop, I don't think you can do it unless there is some structure to your data that can be exploited. The inner loop is very straightforward to vectorize:

for i in xrange(1, d):
    j = np.arange(k-i-1, k+1)
    result[j, i] = ((x - y[j]) / (y[j+i] - y[j]) * result[j, i-1] +
                    (y[j+i+1] - x) / (y[j+i+1] - y[j+1]) * result[j+1, i-1])

There is some marginal improvement you could get by defining a base_j = np.arange(k-d-2, k+1) array outside the loop, and then slicing it inside the loop with something like j = base_j[d-i+1:].

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

Comments

4

You can get rid of the inner loop via clever broadcasting

tmp = ((x - y[k-i+1:k+1])/(y[k+1:k+1+i] - y[j]))*result[k-i+1:k+1,i-1] + ((y[k+2:k+i+2]-x)/(y[k+2:k+i+2] - y[k-i+2:k+2]))*result[k-i+2:k+2,i-1]
result[k-i+1:k+1, i] = tmp

but because you bounds in the inner loop depend on the outer-loop you can not remove it through broadcasting as well.

I used tmp in an overabundance of caution

2 Comments

unclear why you accepted mine, @Jaime's is much more readable.
Yoiu did manage to get it working with slicing instead of fancy indexing, so there's no data copying happening. It should be better performance wise,

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.