Vectorize nested for loop Python

For your question to make sense to me, I think you need to move where list = [] appears. Otherwise you'll never get to even i=0, j=1 until list is full. I can't imagine that it is slow as currently implemented --- list will be full very quickly, and then the for loops should be very fast. Here is what I believe you intended. Please clarify if this is not correct.

for i in range (0,array.shape[0]):
    for j in range (0,array.shape[1]):
         list = []
         while len(list) < 100:
                print "identity", i, j

                #find neighboring entry with greatest value (e.g., assume it is [i-1, j] with value 10)
                list.append((i-1,j))
                i = i-1
                j = j
         #perform operations on list

Let's do some modifications. I'll assume there is a function get_max_nbr(i,j) which returns the coordinates of the maximum neighbor. One of the places your code is slow is that it will call get_max_nbr for the same coordinate many times (at each step in the loop it does it 100 times). The code below uses memoization to get around this (down to 1 time on average). So if this is your bottleneck, this should get you close to 100 times speedup.

maxnbr = {}
for i in range(0,array.shape[0]):
    for j in range (0,array.shape[1]):
        list = []
        current_loc = (i,j)
        while len(list) < 100:
            if current_loc not in maxnbr:  #if this is our first time seeing current_loc
                maxnbr[current_loc] = get_max_nbr(*current_loc) #note func(*(i,j)) becomes func(i,j)
            current_loc = maxnbr[current_loc]
            list.append(current_loc)
        #process list here                 

This doesn't successfully vectorize, but it does create the list (I think) you want, and it should be a significant improvement. It may be that if we knew more about the list processing we might be able to find a better approach, but it's not clear.

Very simply, numpy allows for element-wise operation on its arrays without having to loop over each of its dimensions.

So say you want to apply a simple operator on each element, e.g. scalar multiplication by a number 2, then you can do either of the following:

array*2

or

np.multiply( array,2)

Depending on the nature of stuff you are doing within your loop, you may adapt other techniques to do an elementwise operation using vectorization.

• Your first concern should be to see if you can carry out your calculations using numpy's element-wise operators.
• If that doesn't work, look at the universal functions (ufuncs) built into numpy.

Both of these are coded in compiled C (or Fortran) and are much faster than looping in Python. Additionally, your code will be shorter and easier to understand.

Additional parameters that might improve performance are which compiler was used to compile numpy and which lineair algebra library is used (assuming your code uses linear algebra). E.g. the ATLAS are automatically tuned for the machine they are built on. Intel sells a Fortran compiler and math libraries that are supposed to be very fast on an Intel processor. IIRC, they also parallelize over all available cores.

If your math libraries don't use multiple cores automatically, using the multiprocessing module might be an option. Assuming the problem can be parallelized, this can reduce the runtime (almost) by a factor of 1/N where N is the number of cores. Minus of course the overhead necessary to distribute the problem and gather the results.

Alternatively, for problems that can be parallelized, you could use pyCUDA with numpy if you have a NVidia video card.

If your goal is to find the local maxima in the array, you can use scipy.ndimage.filters.maximum_filter with a 3×3 window and then check for equality:

import numpy
import scipy
import scipy.ndimage

arr = numpy.array([[1, 1, 2, 8],
                   [5, 5, 4, 1],
                   [9, 5, 8, 8],
                   [7, 3, 6, 6]])
maxima = zip(*(scipy.ndimage.filters.maximum_filter(arr, 3) == arr).nonzero())

The speed of this will depend greatly on whether you really need to use only the first 100 and how many maxima there are. If so, breaking out early will probably be faster. Refining your question with the real meat of what you are doing will help us to get a better solution, though.

Adding to the already good answers, here is the commented and fast version to get everything in a list:

import numpy as np
import scipy.ndimage as ndi

#Data generation
data=np.random.randint(100, size=(2000, 2000))
#Maximum extraction using a 3x3 kernel
b=ndi.filters.maximum_filter(data,3) 
#Getting the first 100 entries of b as a 1-D array
max_list=b.flatten()[0:99]

In my test, this code took about 0.2 seconds including the data generation on my Intel i7 CPU and about 3 seconds, when the size of the array is 20k*2k. Timing seems to be no issue here, as I run into memory issues before execution time goes up noticeably.

Nevertheless, you can subdivide the exact same approach into smaller subarrays for larger amounts of data. Keep in mind, that at some point the data handling will take more time than the computation itself.