Optimizing a nested-loop operation

Question

I have this datacube containing data for each pixel of an image (pretty much like hyperspectral imaging). I'm trying to fit a line on each pixel of the image, in an efficient way. Right now, I do it like this:

My datacube is a 6X1024x1024 numpy array, and I have another variable containing the independent variable of my data.

map = np.zeros((1024,1024))
for i in np.mgrid[1:1024]:
    for j in np.mgrid[1:1024]:
        x = independent_variable # This is my independent variable
        y = spec_cube[:,i,j] # The Y data to be fitted is the power at each scale, for a pixel
        index = polyfit(x,y,1) # Outputs the slope and the offset
        map[i,j] = index[0] # The pixel value is the index

I know that nested for loops are pretty much the worst thing to do in general, but I can't think of a better way.

I tried the following but it gives this error: "ValueError: too many values to unpack"

map = np.zeros((1024,1024))
for i,j in map:
    x = independent_variable # This is my independent variable
    y = spec_cube[:,i,j] # The Y data to be fitted is the power at each scale, for a pixel
    index = polyfit(x,y,1) # Outputs the slope and the offset
    map[i,j] = index[0] # The pixel value is the index

This and related questions demonstrate a common problem with python/numpy, when you have tight inner loops that has no equivalent numpy vector operations then you are essentially stuck no matter how hard you try to optimize. In these types of cases one should seriously consider other alternative such as C-extensions or better yet Cython. — Rabih Kodeih, Jan 21 '14 at 12:59

David Wolever · Answer 1 · 2012-06-07T05:06:15.787

5

One method to speed things up: use itertools.product:

for (i, j) in itertools.product(np.mgrid[1:1024], np.mgrid[1:1024]):
    ... stuff ...

The improvement (Python 2.7.1):

In [2]: def multiline():
   ...:     for i in np.mgrid[1:1024]:
   ...:         for j in np.mgrid[1:1024]:
   ...:             pass
   ...:        

In [3]: def single_line():
   ...:     for i, j in product(np.mgrid[1:1024], np.mgrid[1:1024]):
   ...:         pass
   ...:    

In [4]: from itertools import product

In [5]: %timeit multiline()
10 loops, best of 3: 138 ms per loop

In [6]: %timeit single_line()
10 loops, best of 3: 75.6 ms per loop

edited Jun 07 '12 at 05:06

answered Jun 07 '12 at 02:15

David Wolever

148,955
89
346
502

Thanks, it did simplify the code and gave some speed bost. It's not significant enough for my purpose, it takes 412 sec for something I want to execute every 3 seconds... I will have to find a way to avoid this computation. – PhilMacKay Jun 08 '12 at 02:51

score 3 · Accepted Answer · answered Jun 08 '12 at 11:51

Since the operation inside the loop was to find the slope of a line, I went with a less precise method, but using array operations. Basically, to find the slope I did: delta Y / delta X for each adjacent points, then averaged all the slopes.

Turns out it takes a fraction of a second.

Here's the new code:

map = np.zeros((spec_cube.shape[1],spec_cube.shape[2])) # This will be the power index map
x = scale_array
for i in np.mgrid[1:spec_cupe.shape[0]]:
  spec_cube[i-1] = (spec_cube[i]-spec_cube[i-1])/(scale_array[i]-scale_array[i-1])
  map += spec_cube[i-1]
map /= (spec_cube.shape[0]-1)

My script went from 420s to 9.66s!

Congrats on the fix! When you are able, please make sure to mark your answer as 'accepted' so that others will be able to learn from your success. Cheers~ — Andrew Kozak, Jun 08 '12 at 13:53

Optimizing a nested-loop operation

2 Answers2