Jak mogę wektoryzować tę potrójną pętlę przez tablice 2d w numpy?

Question

mogę wyeliminować wszystkie pętle Pythona w tym obliczeń:

result[i,j,k] = (x[i] * y[j] * z[k]).sum() 

gdzie x[i], y[j], z[k] są wektory o długości N i x, y, z ma pierwsze wymiary o długości A, B, C S.t. wyjście ma kształt (A,B,C), a każdy element to suma potrójnego produktu (elementowo).

Mogę go pobrać z 3 do 1 pętli (kod poniżej), ale utknąłem próbując usunąć z ostatniej pętli .

W razie potrzeby mogę wykonać A=B=C (przez niewielką ilość wypełnienia).

# Example with 3 loops, 2 loops, 1 loop (testing omitted) 

N = 100 # more like 100k in real problem 
A = 2 # more like 20 in real problem 
B = 3 # more like 20 in real problem 
C = 4 # more like 20 in real problem 

import numpy 
x = numpy.random.rand(A, N) 
y = numpy.random.rand(B, N) 
z = numpy.random.rand(C, N) 

# outputs of each variant 
result_slow = numpy.empty((A,B,C)) 
result_vec_C = numpy.empty((A,B,C)) 
result_vec_CB = numpy.empty((A,B,C)) 

# 3 nested loops 
for i in range(A): 
    for j in range(B): 
     for k in range(C): 
      result_slow[i,j,k] = (x[i] * y[j] * z[k]).sum() 

# vectorize loop over C (2 nested loops) 
for i in range(A): 
    for j in range(B): 
     result_vec_C[i,j,:] = (x[i] * y[j] * z).sum(axis=1) 

# vectorize one C and B (one loop) 
for i in range(A): 
    result_vec_CB[i,:,:] = numpy.dot(x[i] * y, z.transpose()) 

numpy.testing.assert_almost_equal(result_slow, result_vec_C) 
numpy.testing.assert_almost_equal(result_slow, result_vec_CB) 

Answer 1

Jeśli używasz numpy> 1.6, jest niesamowite np.einsum funkcja:

np.einsum('im,jm,km->ijk',x,y,z) 

co jest równoważne do zapętlonych wersjach. Nie jestem pewien, jak to będzie sprawiedliwe pod względem wydajności, gdy dojdziesz do rozmiaru twoich tablic w rzeczywistym problemie (w rzeczywistości dostaję segfault na mojej maszynie, kiedy przejdę do tych rozmiarów). Drugim rozwiązaniem, które często preferuję dla tego rodzaju problemów, jest ponowne napisanie metody z użyciem cythonu.

Answer 2

Korzystanie z einsum ma wiele sensu w twoim przypadku; ale możesz to zrobić całkiem łatwo ręcznie. Sztuką jest sprawić, że tablice będą nadawać się do siebie nawzajem. Oznacza to zmianę ich kształtu tak, aby każda macierz zmieniała się niezależnie wzdłuż własnej osi. Następnie pomnóż je razem, pozwalając numpy zająć się transmisją; a następnie sumuje się wzdłuż ostatniej (prawej) osi.

>>> x = numpy.arange(2 * 4).reshape(2, 4) 
>>> y = numpy.arange(3 * 4).reshape(3, 4) 
>>> z = numpy.arange(4 * 4).reshape(4, 4) 
>>> (x.reshape(2, 1, 1, 4) * 
... y.reshape(1, 3, 1, 4) * 
... z.reshape(1, 1, 4, 4)).sum(axis=3) 
array([[[ 36, 92, 148, 204], 
     [ 92, 244, 396, 548], 
     [ 148, 396, 644, 892]], 

     [[ 92, 244, 396, 548], 
     [ 244, 748, 1252, 1756], 
     [ 396, 1252, 2108, 2964]]]) 

Można zrobić to nieco bardziej uogólnione za pomocą notacji plasterka wartość newaxis (która jest równa None, więc poniżej będzie działać z None również), a fakt, że sum przyjmuje wartości ujemne osi (z -1 oznaczającym ostatni, -2 oznaczający następny do końca, i tak dalej). W ten sposób nie musisz znać oryginalnego kształtu macierzy; tak długo, jak ich ostatnie osie są kompatybilne, będzie to transmitować pierwsze trzy razem:

>>> (x[:, numpy.newaxis, numpy.newaxis, :] * 
... y[numpy.newaxis, :, numpy.newaxis, :] * 
... z[numpy.newaxis, numpy.newaxis, :, :]).sum(axis=-1) 
array([[[ 36, 92, 148, 204], 
     [ 92, 244, 396, 548], 
     [ 148, 396, 644, 892]], 

     [[ 92, 244, 396, 548], 
     [ 244, 748, 1252, 1756], 
     [ 396, 1252, 2108, 2964]]])