矢量化嵌套for循环Python

Question

我有一个正在迭代的 numpy 数组：

import numpy
import math
array = numpy.array([[1, 1, 2, 8, 2, 2],
               [5, 5, 4, 1, 3, 2],
               [5, 5, 4, 1, 3, 2],
               [5, 5, 4, 1, 3, 2],
               [9, 5, 8, 8, 2, 2],
               [7, 3, 6, 6, 2, 2]])


Pixels = ['U','D','R','L','UL','DL','UR','DR']

for i in range (1,array.shape[0]-1):
    for j in range (1,array.shape[1]-1):


         list = []
         while len(list) < 2:
                iToMakeList = i
                jToMakeList = j

                if iToMakeList > array.shape[0]-1 or iToMakeList < 1 or jToMakeList> array.shape[0]-1 or jToMakeList < 1:

                    break

                PixelCoord = {
            'U' : (iToMakeList-1,jToMakeList),
            'D' : (iToMakeList+1,jToMakeList),
            'R' : (iToMakeList,jToMakeList+1),
            'L' : (iToMakeList,jToMakeList-1),
            'UL' : (iToMakeList-1,jToMakeList-1),
            'DL' : (iToMakeList+1,jToMakeList-1),
            'UR' : (iToMakeList-1,jToMakeList+1),
            'DR' : (iToMakeList+1,jToMakeList+1)
                }
                Value = {
            'U' : array[iToMakeList-1][jToMakeList],
            'D' : array[iToMakeList+1][jToMakeList],
            'R' : array[iToMakeList][jToMakeList+1],
            'L' : array[iToMakeList][jToMakeList-1],
            'UL' : array[iToMakeList-1][jToMakeList-1],
            'DL' : array[iToMakeList+1][jToMakeList-1],
            'UR' : array[iToMakeList-1][jToMakeList+1],
            'DR' : array[iToMakeList+1][jToMakeList+1]
                }


                candidates = []
                for pixel in Pixels:
                    candidates.append((Value[pixel],pixel))

                Lightest = max(candidates)


                list.append(PixelCoord[Lightest[1]])

                iToMakeList = PixelCoord[Lightest[1]][0]
                jToMakeList = PixelCoord[Lightest[1]][1]

我想加快这个过程。它很慢。

假设这个 sn-p 的输出是我的最终目标，我唯一想做的就是加速这个代码。

Answer 1

为了让您的问题对我有意义，我认为您需要移动到出现 list = [] 的位置。否则，在list 已满之前，您将永远无法到达i=0、j=1。我无法想象它现在实现的速度很慢——列表很快就会满，然后 for 循环应该很快。这就是我相信你的意图。如果这不正确，请澄清。

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

让我们做一些修改。我假设有一个函数get_max_nbr(i,j) 返回最大邻居的坐标。您的代码很慢的地方之一是它会多次调用 get_max_nbr 以获得相同的坐标（在循环中的每一步它都会调用 100 次）。下面的代码使用 memoization 来解决这个问题（平均减少 1 次）。因此，如果这是您的瓶颈，这应该可以让您获得接近 100 倍的加速。

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                 

这并没有成功矢量化，但它确实创建了你想要的列表（我认为），它应该是一个显着的改进。如果我们对列表处理有更多的了解，可能会找到更好的方法，但目前还不清楚。

Answer 2

所以这是我的并行方法。首先，我创建一个查找表，其中每个像素都显示最近邻最大值的坐标。 对于我的 intel i7 双核 cpu 上的 100*100 矩阵，代码在大约 2 秒内运行。 到目前为止，代码还没有优化，多处理内部的数据处理有点奇怪，肯定可以变得更容易。让我知道，如果这是，你想要什么。到目前为止，代码仅将数据点的坐标添加到列表中，如果您需要这些值，请在适当的点进行更改或仅解析生成的 lines[] 列表。

import numpy
import multiprocessing as mp
import time
start=time.time()
#Getting the number of CPUs present
num_cpu=mp.cpu_count()
#Creation of random data for testing
data=numpy.random.randint(1,30,size=(200,200))
x,y=data.shape
#Padding is introduced to cope with the border of the dataset.
#Change if you want other behaviour like wrapping, reflection etc.
def pad(data):
    '''Can be made faster, by using numpys pad function
    if present'''
    a=numpy.zeros((x+2,y+2))
    a[1:-1,1:-1]=data
    return a
data=pad(data)
#Kernel to get only the neighbours, change that if you only want diagonals or other shapes.
kernel=numpy.array([[1,1,1],[1,0,1],[1,1,1]])
result_list=[]  
#Setting up functions for Parallel Processing  
def log_result(result): 
    result_list.append(result) 
def max_getter(pixel):
    '''As this function is going to be used in a parallel processing environment,
    the data has to exist globally in order not to have to pickle it in the subprocess'''
    temp=data[pixel[0]-1:pixel[0]+2,pixel[1]-1:pixel[1]+2].copy()*kernel
    #Getting the submatrix without the central pixel
    compare=numpy.max(temp)==temp
    coords=numpy.nonzero(compare)
    if len(coords[0])!=1:
        coords=(coords[0][0],coords[1][0])
    #discards every maximum which is not the first. Change if you want.
    #Converting back to global coordinates
    return (pixel,(pixel[0]+(numpy.asscalar(coords[0])-1),pixel[1]+(numpy.asscalar(coords[1])-1)))
    #This assumes, that the maximum is unique in the subset, if this is not the case adjust here
def parallell_max():
    pool = mp.Pool() 
#You can experiment using more cores if you have hyperthreading and it's not correctly found by cpu_count
    for i in range(1,x+1):

        for j in range(1,y+1):

            pool.apply_async(max_getter, args = ((i,j),),callback=log_result) 
    pool.close()
    pool.join() 


#___________START Parallel Processing________
if __name__ == '__main__':
   # directions={}
    parallell_max()
    directions={}
    for i in result_list:
        directions[i[0]]=i[1]
    #Directions is a dictionary-like lookup-table, where every pixel gives the next pixel in the line
    lines=[]
    #The following code can also be easily parallelized as seen above.
    for i in range(1,x+1):
        for j in range(1,y+1):
            line=[]
            first,second=i,j
            for k in range(100):
                line.append((first,second))
                first,second=directions[(first,second)]
            lines.append(line)
    stop=time.time()
    print stop-start

Answer 3

非常简单，numpy 允许对其数组进行逐元素操作，而无需遍历其每个维度。

假设你想对每个元素应用一个简单的运算符，例如scalar multiplication 通过数字 2，然后您可以执行以下任一操作：

array*2

或

np.multiply( array,2)

根据您在循环中所做的stuff 的性质，您可以采用其他技术来使用vectorization 进行元素操作。

Answer 4

• 您首先应该关注是否可以使用 numpy 的逐元素运算符进行计算。
• 如果这不起作用，请查看 numpy 中内置的 通用函数 (ufunc)。

这两者都是用编译后的 C（或 Fortran）编码的，并且比 Python 中的循环快很多。此外，您的代码将更短且更易于理解。

可能提高性能的其他参数是用于编译 numpy 的编译器和使用的线性代数库（假设您的代码使用线性代数）。例如。 ATLAS 会自动针对构建它们的机器进行调整。英特尔销售的 Fortran 编译器和数学库在英特尔处理器上应该非常快。 IIRC，它们还对所有可用内核进行并行处理。

如果您的数学库不自动使用多核，则可以选择使用multiprocessing 模块。假设问题可以并行化，这可以（几乎）将运行时间减少 1/N 倍，其中 N 是内核数。当然减去分配问题和收集结果所需的开销。

或者，对于可以并行化的问题，如果您有 NVidia 视频卡，您可以使用 pyCUDA 和 numpy。

Answer 5

如果您的目标是在数组中找到局部最大值，您可以使用带有 3×3 窗口的 scipy.ndimage.filters.maximum_filter，然后检查是否相等：

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

这个速度很大程度上取决于您是否真的只需要使用前 100 个以及有多少个最大值。如果是这样，早点爆发可能会更快。不过，用您正在做的事情的真正实质来完善您的问题将有助于我们获得更好的解决方案。

Answer 6

除了已经很好的答案之外，这里是评论和快速版本，可以将所有内容放在列表中：

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]

在我的测试中，这段代码大约需要 0.2 秒，包括我的 Intel i7 CPU 上的数据生成和大约 3 秒，当数组大小为 20k*2k 时。时间在这里似乎没有问题，因为我在执行时间明显增加之前遇到了内存问题。

尽管如此，您可以将完全相同的方法细分为更小的子数组以获取大量数据。请记住，在某些时候，数据处理将花费比计算本身更多的时间。