Làm cách nào tôi có thể tối ưu hóa phép tính trên hàm này theo dạng gumpy?

Question

Tôi muốn thực hiện các vấn đề sau đây trong numpy và đây là mã của tôi.

Tôi đã thử mã NumPy sau cho vấn đề này với một vòng lặp for. Tôi tự hỏi liệu có cách nào hiệu quả hơn để thực hiện phép tính này không? Tôi thực sự đánh giá cao điều đó!

k, d = X.shape 
m = Y.shape[0] 

c1 = 2.0*sigma**2 
c2 = 0.5*np.log(np.pi*c1) 
c3 = np.log(1.0/k) 

L_B = np.zeros((m,)) 
for i in xrange(m): 
    if i % 100 == 0: 
     print i 
    L_B[i] = np.log(np.sum(np.exp(np.sum(-np.divide(
       np.power(X-Y[i,:],2), c1)-c2,1)+c3))) 

print np.mean(L_B) 

Tôi đã nghĩ đến việc np.expand_dims(X, 2).repeat(Y.shape[0], 2)-Y bằng cách tạo một tensor 3D do đó việc tính toán sau đây có thể được thực hiện bằng cách phát sóng, nhưng điều đó sẽ lãng phí rất nhiều bộ nhớ khi m là lớn.

Tôi cũng tin rằng np.einsum() không sử dụng gì ngoài vòng lặp vì vậy có thể không hiệu quả, hãy sửa tôi nếu tôi sai.

Bạn nghĩ gì?

Answer 1

Tối ưu hóa Stage # 1

cấp độ đầu tiên của tôi về tối ưu hóa sử dụng một dịch trực tiếp của mã điên rồ đến một broadcasting dựa một khi giới thiệu một trục mới và như vậy người ta không nên nhớ hiệu quả, như được liệt kê dưới đây -

p1 = (-((X[:,None] - Y)**2)/c1)-c2 
p11 = p1.sum(2) 
p2 = np.exp(p11+c3) 
out = np.log(p2.sum(0)).mean() 

Tối ưu hóa Stage # 2

Mang trong vài tối ưu hóa lưu giữ trong tâm trí thứ tại chúng tôi có ý định tách ra hoạt động trên các hằng số, tôi đã kết thúc với những điều sau -

c10 = -c1 
c20 = X.shape[1]*c2 

subs = (X[:,None] - Y)**2 
p00 = subs.sum(2) 
p10 = p00/c10 
p11 = p10-c20 
p2 = np.exp(p11+c3) 
out = np.log(p2.sum(0)).mean() 

Tối ưu hóa Stage # 3

Đi xa hơn với nó và và nhìn thấy nơi các hoạt động có thể là tối ưu hóa, tôi đã sử dụng Scipy's cdist để thay thế công việc nặng nhọc của bình phương và sum-reduction. Điều này cần được bộ nhớ khá hiệu quả và đã cho chúng tôi thực hiện cuối cùng, như hình dưới đây -

from scipy.spatial.distance import cdist 

# Setup constants 
c10 = -c1 
c20 = X.shape[1]*c2 
c30 = c20-c3 
c40 = np.exp(c30) 
c50 = np.log(c40) 

# Get stagewise operations corresponding to loopy ones 
p1 = cdist(X, Y, 'sqeuclidean') 
p2 = np.exp(p1/c10).sum(0) 
out = np.log(p2).mean() - c50 

---

Runtime kiểm tra

phương pháp tiếp cận -

def loopy_app(X, Y, sigma): 
    k, d = X.shape 
    m = Y.shape[0] 

    c1 = 2.0*sigma**2 
    c2 = 0.5*np.log(np.pi*c1) 
    c3 = np.log(1.0/k) 

    L_B = np.zeros((m,)) 
    for i in xrange(m): 
     L_B[i] = np.log(np.sum(np.exp(np.sum(-np.divide(
        np.power(X-Y[i,:],2), c1)-c2,1)+c3))) 

    return np.mean(L_B) 

def vectorized_app(X, Y, sigma): 
    # Setup constants 
    k, d = D_A.shape 
    c1 = 2.0*sigma**2 
    c2 = 0.5*np.log(np.pi*c1) 
    c3 = np.log(1.0/k) 

    c10 = -c1 
    c20 = X.shape[1]*c2 
    c30 = c20-c3 
    c40 = np.exp(c30) 
    c50 = np.log(c40) 

    # Get stagewise operations corresponding to loopy ones 
    p1 = cdist(X, Y, 'sqeuclidean') 
    p2 = np.exp(p1/c10).sum(0) 
    out = np.log(p2).mean() - c50 
    return out 

Thời gian và xác minh -

In [294]: # Setup inputs with m(=D_B.shape[0]) being a large number 
    ...: X = np.random.randint(0,9,(100,10)) 
    ...: Y = np.random.randint(0,9,(10000,10)) 
    ...: sigma = 2.34 
    ...: 

In [295]: np.allclose(loopy_app(X, Y, sigma),vectorized_app(X, Y, sigma)) 
Out[295]: True 

In [296]: %timeit loopy_app(X, Y, sigma) 
1 loops, best of 3: 225 ms per loop 

In [297]: %timeit vectorized_app(X, Y, sigma) 
10 loops, best of 3: 23.6 ms per loop 

In [298]: # Setup inputs with m(=Y.shape[0]) being a much large number 
    ...: X = np.random.randint(0,9,(100,10)) 
    ...: Y = np.random.randint(0,9,(100000,10)) 
    ...: sigma = 2.34 
    ...: 

In [299]: np.allclose(loopy_app(X, Y, sigma),vectorized_app(X, Y, sigma)) 
Out[299]: True 

In [300]: %timeit loopy_app(X, Y, sigma) 
1 loops, best of 3: 2.27 s per loop 

In [301]: %timeit vectorized_app(X, Y, sigma) 
1 loops, best of 3: 243 ms per loop 

Khoảng 10x tăng tốc ở đó!