Как векторизовать массивы 3D Numpy

Question

У меня есть массив 3D numpy, такой как a = np.zeros((100,100, 20)). Я хочу выполнить операцию над каждой позицией x,y, которая включает в себя все элементы над осью z, и результат сохраняется в массиве, таком как b = np.zeros((100,100)), на той же соответствующей позиции x,y.

Теперь я делаю это с помощью для цикла:

d_n = np.array([...]) # a parameter with the same shape as b 
for (x,y), v in np.ndenumerate(b): 
    C = a[x,y,:] 

    ### calculate some_value using C 
    minv = sys.maxint 
    depth = -1 
    C = a[x,y,:] 
    for d in range(len(C)): 
     e = 2.5 * float(math.pow(d_n[x,y] - d, 2)) + C[d] * 0.05 
     if e < minv: 
      minv = e 
      depth = d 

    some_value = depth 
    if depth == -1: 
     some_value = len(C) - 1 
    ### 

    b[x,y] = some_value 

Сейчас проблема заключается в том, что эта операция гораздо медленнее, чем другие выполнявший вещий путь, например c = b * b (я на самом деле профилированный эта функцию, и это примерно на 2 порядка медленнее, чем другие, использующих NumPy встроенных функций и векторизованных функции, по сравнению с аналогичным числом элементов)

Как я могу улучшить производительность такого рода функции, отображающие 3D-массив для двумерного?

Answer 1

что обычно делается в 3D-изображения, чтобы поменять местами оси Z к первому индексу:

>>> a = a.transpose((2,0,1)) 
>>> a.shape 
(20, 100, 100) 

И теперь вы можете легко перебирать по оси Z:

>>> for slice in a: 
     do something 

slice здесь будет каждый из ваших 100x100 фракций 3D-матрицы. Кроме того, благодаря транспозиции вы можете получить доступ к каждому из 2D-фрагментов напрямую, индексируя первую ось. Например, a[10] предоставит вам 11-й 2D 100x100 фрагмент.

Бонус:. Если вы храните данные contiguosly без транспозиции (или преобразования в непрерывный массив, используя a = np.ascontiguousarray(a.transpose((2,0,1))) доступ к вам 2D ломтиков будет быстрее, так как они отображаются contiguosly в памяти

Answer 2

Очевидно, что вы хотите, чтобы избавиться от явного for цикла, но я думаю, возможно ли это зависит от того, что расчет вы делаете с C. В качестве простого примера,

a = np.zeros((100,100, 20)) 
a[:,:] = np.linspace(1,20,20) # example data: 1,2,3,.., 20 as "z" for every "x","y" 

b = np.sum(a[:,:]**2, axis=2) 

заполнит 100 по 100 массива b с суммой квадратов «г» значений a, то есть 1 + 4 + 9 + ... + 400 = 2870.

Answer 3

Если ваш внутренний расчет достаточно сложна, и не поддаются векторизации, то ваша структура итерация хорошо, и не вносит значительного вклада в момент расчета

for (x,y), v in np.ndenumerate(b): 
    C = a[x,y,:] 
    ... 
    for d in range(len(C)): 
     ... # complex, not vectorizable calc 
    ... 
    b[x,y] = some_value 

Там не кажется, специальная структура, в 1-й 2-х измерениях, так что вы может так же хорошо воспринимать это как 2D-отображение на 1D, например карту ping a (N,20) массив на массив (N,). Это ничего не ускоряет, но может помочь выделить существенную структуру проблемы.

Один шаг - сосредоточиться на ускорении этого C до some_value расчет. Существуют такие функции, как cumsum и cumprod, которые помогут вам выполнять последовательные вычисления по вектору. cython - также хороший инструмент.

Другой подход заключается в том, чтобы увидеть, можете ли вы выполнить этот внутренний расчет по значениям N все сразу. Другими словами, если вы должны итерации, лучше сделать это в самом маленьком измерении.

В некотором смысле это не-ответ. Но без полного понимания того, как вы получаете some_value от C и d_n Я не думаю, что мы можем сделать больше.

---

Похоже e могут быть вычислены для всех точек сразу:

e = 2.5 * float(math.pow(d_n[x,y] - d, 2)) + C[d] * 0.05 

E = 2.5 * (d_n[...,None] - np.arange(a.shape[-1]))**2 + a * 0.05 # (100,100,20) 

E.min(axis=-1) # smallest value along the last dimension 
E.argmin(axis=-1) # index of where that min occurs 

На первый взгляд это выглядит как этот E.argmin это значение b, что вы хотите (переделаны для некоторых граничных условий при необходимости).

Я не реалистические a и d_n массивов, но с простым проверочными, это E.argmin(-1) соответствует вашему b с 66x SpeedUp.

Answer 4

Как я могу улучшить производительность таких функций, сопоставляя 3D-массив с двумерным?

Многие функции в Numpy являются «сокращение» функции ^*, например sum, any, std и т.д. Если вы поставлять axis аргумент, кроме None такой функции, что позволит уменьшить размерность массив по этой оси. Для вашего кода вы можете использовать функцию argmin, если вы сначала вычислить e в векторизованном образе:

d = np.arange(a.shape[2]) 
e = 2.5 * (d_n[...,None] - d)**2 + a*0.05 
b = np.argmin(e, axis=2) 

Индексация с [...,None] используются заниматься broadcasting. Значения в e имеют значения с плавающей точкой, так что это немного странно сравнивать с sys.maxint но там вы идете:

I, J = np.indices(b.shape) 
b[e[I,J,b] >= sys.maxint] = a.shape[2] - 1 

_{* Strickly говоря функция снижения имеет вид reduce(operator, sequence) так что технически не std и argmin}