2015-10-08 25 views
9

Tôi nhận được một số kết quả kiểm tra hiệu quả mà tôi không thể giải thích.Tại sao B = numpy.dot (A, x) lặp đi lặp lại chậm hơn nhiều khi thực hiện B [i,:,:] = numpy.dot (A [i,:,:], x))?

Tôi muốn lắp ráp ma trận B có mục nhập thứ i B [i,:,:] = A [i,:,:]. Chấm (x), trong đó mỗi A [i,:,:] là một ma trận 2D, và như vậy là x.

Tôi có thể thực hiện ba cách này, để kiểm tra hiệu suất tôi tạo ngẫu nhiên (numpy.random.randn) ma trận A = (10,1000,1000), x = (1000,1200). Tôi có được kết quả thời gian sau:

(1) đơn đa chiều dot sản phẩm

B = A.dot(x) 

total time: 102.361 s 

(2) vòng lặp thông qua i và thực hiện các sản phẩm chấm 2D

# initialize B = np.zeros([dim1, dim2, dim3]) 
    for i in range(A.shape[0]): 
     B[i,:,:] = A[i,:,:].dot(x) 

total time: 0.826 s 

(3) numpy.einsum

B3 = np.einsum("ijk, kl -> ijl", A, x) 

total time: 8.289 s 

Vì vậy, tùy chọn (2) là nhanh nhất cho đến nay. Nhưng, chỉ xem xét (1) và (2), tôi không thấy sự khác biệt lớn giữa chúng. Làm thế nào có thể lặp đi lặp lại và làm các sản phẩm chấm 2D nhanh hơn ~ 124 lần? Cả hai đều sử dụng numpy.dot. Bất kỳ thông tin chi tiết nào?

tôi bao gồm các mã được sử dụng cho các kết quả trên chỉ dưới đây:

import numpy as np 
import numpy.random as npr 
import time 

dim1, dim2, dim3 = 10, 1000, 1200 
A = npr.randn(dim1, dim2, dim2) 
x = npr.randn(dim2, dim3) 

# consider three ways of assembling the same matrix B: B1, B2, B3 

t = time.time() 
B1 = np.dot(A,x) 
td1 = time.time() - t 
print "a single dot product of A [shape = (%d, %d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td1) 


B2 = np.zeros([A.shape[0], x.shape[0], x.shape[1]]) 
t = time.time() 
for i in range(A.shape[0]): 
    B2[i,:,:] = np.dot(A[i,:,:], x) 
td2 = time.time() - t 
print "taking %d dot products of 2D dot products A[i,:,:] [shape = (%d, %d)] with x [shape = (%d, %d)] completes in %.3f s" \ 
    % (A.shape[0], A.shape[1], A.shape[2], x.shape[0], x.shape[1], td2) 

t = time.time() 
B3 = np.einsum("ijk, kl -> ijl", A, x) 
td3 = time.time() - t 
print "using np.einsum, it completes in %.3f s" % td3 

Trả lời

3

Với dims nhỏ 10,100,200, tôi nhận được một thứ hạng tương tự

In [355]: %%timeit 
    .....: B=np.zeros((N,M,L)) 
    .....: for i in range(N): 
       B[i,:,:]=np.dot(A[i,:,:],x) 
    .....: 
10 loops, best of 3: 22.5 ms per loop 
In [356]: timeit np.dot(A,x) 
10 loops, best of 3: 44.2 ms per loop 
In [357]: timeit np.einsum('ijk,km->ijm',A,x) 
10 loops, best of 3: 29 ms per loop 

In [367]: timeit np.dot(A.reshape(-1,M),x).reshape(N,M,L) 
10 loops, best of 3: 22.1 ms per loop 

In [375]: timeit np.tensordot(A,x,(2,0)) 
10 loops, best of 3: 22.2 ms per loop 

các itererative là nhanh hơn, mặc dù không nhiều như trong trường hợp của bạn.

Điều này có thể đúng miễn là thứ nguyên lặp đó nhỏ hơn so với các thứ nguyên khác. Trong trường hợp đó, chi phí lặp lại (các cuộc gọi hàm vv) là nhỏ so với thời gian tính toán. Và làm tất cả các giá trị cùng một lúc sẽ sử dụng nhiều bộ nhớ hơn.

Tôi đã thử một biến thể dot nơi tôi định hình lại A thành 2d, nghĩ rằng dot thực hiện loại định hình lại nội bộ. Tôi ngạc nhiên rằng nó thực sự là nhanh nhất. tensordot có thể thực hiện cùng một định dạng lại (mã đó nếu Python có thể đọc được).


einsum thiết lập một 'tổng hợp các sản phẩm' lặp liên quan đến 4 biến, i,j,k,m - đó là dim1*dim2*dim2*dim3 bước với mức độ C nditer. Vì vậy, càng có nhiều chỉ số bạn có không gian lặp lớn hơn.

+0

'dấu chấm' thực hiện nhiều điều dưới sự che chở, rõ ràng là' np.dot (A, x) 'không gọi BLAS và bằng cách nào đó là mặc định đối với thường trình GEMM bên trong của numpy. Ví dụ reshape của bạn là bỏ qua cơ chế looping và đi thẳng đến một cuộc gọi GEMM 2D thông thường mà không cần sao chép bất kỳ dữ liệu nào, nó sẽ luôn là giải pháp nhanh nhất cho các vấn đề có kích thước hợp lý cho một BLAS tốt. Trên máy tính xách tay của tôi với MKL của nó ~ 50 lần nhanh hơn einsum cho các vấn đề kích thước ban đầu. – Daniel

+0

'tensordot' đang thực hiện cùng một định dạng lại. – hpaulj

+0

Vâng tensordot buộc một bản sao của dữ liệu nội bộ, tôi sẽ không khuyên bạn nên tùy chọn này. – Daniel

2

Tôi không quá quen thuộc với C-API NumPy, và các numpy.dot là một trong những chức năng được xây dựng trong đó được sử dụng để dưới _dotblas ở trước phiên bản.

Tuy nhiên, đây là suy nghĩ của tôi.

1)numpy.dot có các đường dẫn khác nhau cho mảng 2 chiều và mảng n chiều. Từ 's online documentationnumpy.dot:

Đối với mảng 2-D nó tương đương với phép nhân ma trận, và đối với mảng 1 chiều đến sản phẩm bên trong của vectơ (không chia phức tạp). Đối với các thứ nguyên N, nó là một sản phẩm tổng hợp trên trục cuối cùng của một và số thứ hai đến cuối cùng của b

dấu chấm (a, b) [i, j, k, m] = sum (a [i, j ,:] * b [k,:, m])

vì vậy, đối với mảng 2-D bạn luôn được đảm bảo để có một cuộc gọi đến BLAS của dgemm, tuy nhiên đối với mảng NĐ numPy có thể chọn nhân trục cho mảng đó có thể không tương ứng với trục thay đổi nhanh nhất (như bạn có thể thấy từ đoạn trích tôi đã đăng), và kết quả là toàn bộ sức mạnh của dgemm có thể bị bỏ qua.

2) mảng A của bạn quá lớn để được tải vào bộ đệm CPU. Trong ví dụ của bạn, bạn sử dụng A với kích thước (10,1000,1000) mang đến cho

In [1]: A.nbytes 
80000000 
In [2]: 80000000/1024 
78125 

Đó là gần như 80MB, lớn hơn nhiều so với kích thước bộ nhớ cache. Vì vậy, một lần nữa bạn mất hầu hết quyền lực của dgemm ngay tại đó.

3) Bạn cũng đang định thời gian các chức năng hơi không chính xác. Hàm time trong Python được biết là không chính xác. Sử dụng timeit để thay thế.

Vì vậy, có tất cả các điểm trên trong tâm trí, chúng ta hãy cố gắng thử nghiệm với các mảng có thể được nạp vào bộ nhớ cache

dim1, dim2, dim3 = 20, 20, 20 
A = np.random.rand(dim1, dim2, dim2) 
x = np.random.rand(dim2, dim3) 

def for_dot1(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[i,:,:], x) 

def for_dot2(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,i,:], x)  

def for_dot3(A,x): 
    for i in range(A.shape[0]): 
     np.dot(A[:,:,i], x) 

và đây là timings mà tôi nhận được (sử dụng numpy 1.9.2 xây dựng chống OpenBLAS 0.2.14):

In [3]: %timeit np.dot(A,x) 
10000 loops, best of 3: 174 µs per loop 
In [4]: %timeit np.einsum("ijk, kl -> ijl", A, x) 
10000 loops, best of 3: 108 µs per loop 
In [5]: %timeit np.einsum("ijk, lk -> ijl", A, x) 
10000 loops, best of 3: 97.1 µs per loop 
In [6]: %timeit np.einsum("ikj, kl -> ijl", A, x) 
1000 loops, best of 3: 238 µs per loop 
In [7]: %timeit np.einsum("kij, kl -> ijl", A, x) 
10000 loops, best of 3: 113 µs per loop 
In [8]: %timeit for_dot1(A,x) 
10000 loops, best of 3: 101 µs per loop 
In [9]: %timeit for_dot2(A,x) 
10000 loops, best of 3: 131 µs per loop 
In [10]: %timeit for_dot3(A,x) 
10000 loops, best of 3: 133 µs per loop 

Lưu ý rằng vẫn còn chênh lệch thời gian, nhưng không theo thứ tự độ lớn. Cũng lưu ý tầm quan trọng của choosing the axis of multiplication. Bây giờ có lẽ, một nhà phát triển vất vả có thể làm sáng tỏ những gì mà numpy.dot thực sự thực hiện dưới sự che chở cho các mảng N-D.

+0

Đây là 1). 'time.time()' là hợp lệ miễn là bạn không xử lý các hàm có độ dài micro/nano giây. Trong ví dụ của bạn, bạn hiển thị đối số trục chỉ là hệ số 2/3 trong khi các vấn đề thời gian khác nhau theo hệ số 1000. Ngoài ra, hãy lưu ý rằng đối với nhà cung cấp BLAS GEMM (N^3 tính toán, N^2 dữ liệu), bộ đệm nên có phương sai tương đối nhỏ. – Daniel

+0

Vâng, đối với bộ nhớ đệm BLAS của nhà cung cấp có ít ảnh hưởng. Cảm ơn, rõ ràng bây giờ, lý do thực tế cho vấn đề thời gian là do Numpy gọi gemm nội bộ của nó chứ không phải là nhà cung cấp BLAS. – romeric

+0

Tôi nghĩ [this] (https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/multiarraymodule.c#L986) là dòng có liên quan trong nguồn C. Nó được gọi thông qua [chức năng này] (https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/methods.c#L1991). –

Các vấn đề liên quan