Có cách nào để làm cho numpy.argmin() nhanh như min() không?

Question

Tôi đang cố gắng tìm các chỉ số mảng tối thiểu dọc theo một chiều của một mảng có khối lượng 2D rất lớn. Tôi thấy rằng điều này là rất chậm (đã cố gắng tăng tốc nó với nút cổ chai, mà chỉ là một cải tiến tối thiểu). Tuy nhiên, lấy tối thiểu thẳng dường như là một thứ tự cường độ nhanh:

import numpy as np 
import time 

randvals = np.random.rand(3000,160000) 
start = time.time() 
minval = randvals.min(axis=0) 
print "Took {0:.2f} seconds to compute min".format(time.time()-start) 
start = time.time() 
minindex = np.argmin(randvals,axis=0) 
print "Took {0:.2f} seconds to compute argmin".format(time.time()-start) 

Trên máy tính của tôi kết quả đầu ra này:

Took 0.83 seconds to compute min 
Took 9.58 seconds to compute argmin 

Có bất kỳ lý do tại sao argmin rất chậm hơn? Có cách nào để tăng tốc độ để so sánh với min?

Answer 1

In [1]: import numpy as np 

In [2]: a = np.random.rand(3000, 16000) 

In [3]: %timeit a.min(axis=0) 
1 loops, best of 3: 421 ms per loop 

In [4]: %timeit a.argmin(axis=0) 
1 loops, best of 3: 1.95 s per loop 

In [5]: %timeit a.min(axis=1) 
1 loops, best of 3: 302 ms per loop 

In [6]: %timeit a.argmin(axis=1) 
1 loops, best of 3: 303 ms per loop 

In [7]: %timeit a.T.argmin(axis=1) 
1 loops, best of 3: 1.78 s per loop 

In [8]: %timeit np.asfortranarray(a).argmin(axis=0) 
1 loops, best of 3: 1.97 s per loop 

In [9]: b = np.asfortranarray(a) 

In [10]: %timeit b.argmin(axis=0) 
1 loops, best of 3: 329 ms per loop 

lẽ min là đủ thông minh để làm tuần tự công việc của mình trong mảng (do đó với địa phương cache), và argmin là nhảy quanh mảng (gây ra rất nhiều bỏ lỡ bộ nhớ cache)?

Dù sao, nếu bạn sẵn sàng giữ randvals làm mảng được đặt hàng trước Fortran ngay từ đầu, nó sẽ nhanh hơn, mặc dù việc sao chép vào trong Fortran không có tác dụng.

Answer 2

Tôi chỉ cần lấy một cái nhìn tại the source code, và trong khi tôi không hoàn toàn hiểu tại sao mọi thứ đang được thực hiện theo cách họ đang có, đây là những gì sẽ xảy ra:

1. np.min về cơ bản là một cuộc gọi đến np.minimum.reduce.

2. np.argmin di chuyển đầu tiên trên trục bạn muốn hoạt động trên đến hết các tuple hình dạng, sau đó làm cho nó một mảng liền kề, trong đó tất nhiên gây ra một bản sao của mảng đầy đủ trừ khi trục là chiếc cuối cùng để bắt đầu với .

Từ một bản sao đang được thực hiện, bạn có thể sáng tạo và cố gắng nhanh chóng mảng rẻ:

a = np.random.rand(1000, 2000) 

def fast_argmin_axis_0(a): 
    matches = np.nonzero((a == np.min(a, axis=0)).ravel())[0] 
    rows, cols = np.unravel_index(matches, a.shape) 
    argmin_array = np.empty(a.shape[1], dtype=np.intp) 
    argmin_array[cols] = rows 
    return argmin_array 

In [8]: np.argmin(a, axis=0) 
Out[8]: array([230, 532, 815, ..., 670, 702, 989], dtype=int64) 

In [9]: fast_argmin_axis_0(a) 
Out[9]: array([230, 532, 815, ..., 670, 702, 989], dtype=int64) 

In [10]: %timeit np.argmin(a, axis=0) 
10 loops, best of 3: 27.3 ms per loop 

In [11]: %timeit fast_argmin_axis_0(a) 
100 loops, best of 3: 15 ms per loop 

tôi sẽ không đi xa như kêu gọi việc thực hiện hiện một lỗi, vì có thể là lý do chính đáng để làm những gì nó làm theo cách của nó, nhưng loại thủ thuật này có thể tăng tốc độ chức năng được tối ưu hóa cao, mạnh mẽ gợi ý rằng mọi thứ có thể được thực hiện tốt hơn.

Answer 3

Tôi vừa mới gặp vấn đề tương tự và thấy sự khác biệt lớn về hiệu suất khi trục 0 được chọn để tìm giá trị tối thiểu. Như đã đề xuất, vấn đề dường như có liên quan đến việc sao chép nội bộ mảng.

Tôi nghĩ ra một cách giải quyết khá đơn giản bằng cách sử dụng Cython đồng thời thiết lập các giá trị tối thiểu và vị trí của chúng trong một mảng numpy 2D dọc theo một trục nhất định.Lưu ý rằng đối axis = 0, thuật toán hoạt động trên một loạt các cột (số lượng tối đa theo quy định của hằng số blocksize - đây thiết lập tương đương với 8 Kbyte dữ liệu) đồng thời để tận dụng tốt các bộ nhớ cache:

%%cython -c=-O2 -c=-march=native 

import numpy as np 
cimport numpy as np 
cimport cython 
from libc.stdint cimport uint32_t 

@cython.boundscheck(False) 
@cython.wraparound(False) 
cdef void _minargmin_2d_64_axis_0(uint32_t[:] minloc, double[:] minimum, double[:, :] arr) nogil: 
    """ 
    Find the minimum and argmin for a 2D array of 64-bit floats along axis 0 
    Parameters: 
    ----------- 
    arr: numpy_array, dtype=np.float64, shape=(m, n) 
     The array for which the minima should be computed. 
    minloc: numpy_array, dtype=np.uint32, shape=(n) 
     Stores the rows where the minima occur for each column. 
    minimum: numpy_array, dtype=np.float64, shape=(n) 
     The minima along each column. 
    """ 

    # Columns of the matrix are accessed in blocks to increase cache performance. 
    # Specify the number of columns here: 

    DEF blocksize = 1024 

    cdef int i, j, k 
    cdef double minim[blocksize] 
    cdef uint32_t minimloc[blocksize] 

    # Read blocks of data to make good use of the cache 

    cdef int blocks 
    blocks = arr.shape[1]/blocksize 
    remains = arr.shape[1] % blocksize 

    for i in xrange(0, blocksize * blocks, blocksize): 

     for k in xrange(blocksize): 
      minim[k] = arr[0, i + k] 
      minimloc[k] = 0 

     for j in xrange(1, arr.shape[0]): 

      for k in xrange(blocksize): 

       if arr[j, i + k] < minim[k]: 
        minim[k] = arr[j, i + k] 
        minimloc[k] = j 

     for k in xrange(blocksize): 
      minimum[i + k] = minim[k] 
      minloc[i + k] = minimloc[k] 

    # Work on the final 'incomplete' block 

    i = blocksize * blocks 

    for k in xrange(remains): 
     minim[k] = arr[0, i + k] 
     minimloc[k] = 0 

    for j in xrange(1, arr.shape[0]): 

     for k in xrange(remains): 

      if arr[j, i + k] < minim[k]: 
       minim[k] = arr[j, i + k] 
       minimloc[k] = j 

    for k in xrange(remains): 
     minimum[i + k] = minim[k] 
     minloc[i + k] = minimloc[k] 

    # Done! 


@cython.boundscheck(False) 
@cython.wraparound(False) 
cdef void _minargmin_2d_64_axis_1(uint32_t[:] minloc, double[:] minimum, double[:, :] arr) nogil: 
    """ 
    Find the minimum and argmin for a 2D array of 64-bit floats along axis 1 
    Parameters: 
    ----------- 
    arr: numpy_array, dtype=np.float64, shape=(m, n) 
     The array for which the minima should be computed. 
    minloc: numpy_array, dtype=np.uint32, shape=(m) 
     Stores the rows where the minima occur for each row. 
    minimum: numpy_array, dtype=np.float64, shape=(m) 
     The minima along each row. 
    """ 
    cdef int i 
    cdef int j 
    cdef double minim 
    cdef uint32_t minimloc 

    for i in xrange(arr.shape[0]): 
     minim = arr[i, 0] 
     minimloc = 0 

     for j in xrange(1, arr.shape[1]): 
      if arr[i, j] < minim: 
       minim = arr[i, j] 
       minimloc = j 

     minimum[i] = minim 
     minloc[i] = minimloc 


@cython.boundscheck(False) 
@cython.wraparound(False) 
cdef void _minargmin_2d_32_axis_0(uint32_t[:] minloc, float[:] minimum, float[:, :] arr) nogil: 
    """ 
    Find the minimum and argmin for a 2D array of 32-bit floats along axis 0 
    Parameters: 
    ----------- 
    arr: numpy_array, dtype=np.float32, shape=(m, n) 
     The array for which the minima should be computed. 
    minloc: numpy_array, dtype=np.uint32, shape=(n) 
     Stores the rows where the minima occur for each column. 
    minimum: numpy_array, dtype=np.float32, shape=(n) 
     The minima along each column. 
    """ 

    # Columns of the matrix are accessed in blocks to increase cache performance. 
    # Specify the number of columns here: 

    DEF blocksize = 2048 

    cdef int i, j, k 
    cdef float minim[blocksize] 
    cdef uint32_t minimloc[blocksize] 

    # Read blocks of data to make good use of the cache 

    cdef int blocks 
    blocks = arr.shape[1]/blocksize 
    remains = arr.shape[1] % blocksize 

    for i in xrange(0, blocksize * blocks, blocksize): 

     for k in xrange(blocksize): 
      minim[k] = arr[0, i + k] 
      minimloc[k] = 0 

     for j in xrange(1, arr.shape[0]): 

      for k in xrange(blocksize): 

       if arr[j, i + k] < minim[k]: 
        minim[k] = arr[j, i + k] 
        minimloc[k] = j 

     for k in xrange(blocksize): 
      minimum[i + k] = minim[k] 
      minloc[i + k] = minimloc[k] 

    # Work on the final 'incomplete' block 

    i = blocksize * blocks 

    for k in xrange(remains): 
     minim[k] = arr[0, i + k] 
     minimloc[k] = 0 

    for j in xrange(1, arr.shape[0]): 

     for k in xrange(remains): 

      if arr[j, i + k] < minim[k]: 
       minim[k] = arr[j, i + k] 
       minimloc[k] = j 

    for k in xrange(remains): 
     minimum[i + k] = minim[k] 
     minloc[i + k] = minimloc[k] 

    # Done! 

@cython.boundscheck(False) 
@cython.wraparound(False) 
cdef void _minargmin_2d_32_axis_1(uint32_t[:] minloc, float[:] minimum, float[:, :] arr) nogil: 
    """ 
    Find the minimum and argmin for a 2D array of 32-bit floats along axis 1 
    Parameters: 
    ----------- 
    arr: numpy_array, dtype=np.float32, shape=(m, n) 
     The array for which the minima should be computed. 
    minloc: numpy_array, dtype=np.uint32, shape=(m) 
     Stores the rows where the minima occur for each row. 
    minimum: numpy_array, dtype=np.float32, shape=(m) 
     The minima along each row. 
    """ 
    cdef int i 
    cdef int j 
    cdef float minim 
    cdef uint32_t minimloc 

    for i in xrange(arr.shape[0]): 
     minim = arr[i, 0] 
     minimloc = 0 

     for j in xrange(1, arr.shape[1]): 
      if arr[i, j] < minim: 
       minim = arr[i, j] 
       minimloc = j 

     minimum[i] = minim 
     minloc[i] = minimloc 

def Min_Argmin(array_2d, axis = 1): 
    """ 
    Find the minima and corresponding locations (argmin) of a two-dimensional 
    numpy array of floats along a given axis simultaneously, and returns them 
    as a tuple of arrays (min_2d, argmin_2d). 

    (Note: float16 arrays will be internally transformed to float32 arrays.) 
    ---------- 
    array_2d: numpy_array, dtype=np.float32 or np.float64, shape=(m, n) 
     The array for which the minima should be computed. 
    axis : int, 0 or 1 (default) : 
     The axis along which minima are computed. 
    min_2d: numpy_array, dtype=np.uint8, shape=(m) or shape=(n): 
     The array where the minima along the given axis are stored. 
    argmin_2d: 
     The array storing the row/column where the minimum occurs. 
    """ 

    # Sanity checks: 
    if len(array_2d.shape) != 2: 
     raise IndexError('Min_Argmin: Number of dimensions of array must be 2') 

    if not (axis == 0 or axis == 1): 
     raise ValueError('Min_Argmin: Invalid axis specified') 

    arr_type = array_2d.dtype 

    if not arr_type in ('float16', 'float32', 'float64'): 
     raise ValueError('Min_Argmin: Not a float array') 

    # Transform float16 arrays 
    if arr_type == 'float16': 
     array_2d = array_2d.astype(np.float32) 

    # Run analysis 

    if arr_type == 'float64': 

     # Double accuracy 

     min_array = np.zeros(array_2d.shape[1 - axis], dtype = np.float64) 
     argmin_array = np.zeros(array_2d.shape[1 - axis], dtype = np.uint32) 

     if (axis == 0): 
      _minargmin_2d_64_axis_0(argmin_array, min_array, array_2d) 

     else: 
      _minargmin_2d_64_axis_1(argmin_array, min_array, array_2d) 

    else: 

     # Single accuracy 

     min_array = np.zeros(array_2d.shape[1 - axis], dtype = np.float32) 
     argmin_array = np.zeros(array_2d.shape[1 - axis], dtype = np.uint32) 

     if (axis == 0): 
      _minargmin_2d_32_axis_0(argmin_array, min_array, array_2d) 

     else: 
      _minargmin_2d_32_axis_1(argmin_array, min_array, array_2d) 

     # Transform back to float16 type as necessary 

     if arr_type == 'float16': 
      min_array = min_array.astype(np.float16) 


    # Return results 
    return min_array, argmin_array 

Mã này có thể được đặt và biên soạn trong một tế bào máy tính xách tay IPython sau khi tải hỗ trợ Cython:

%load_ext Cython 

và sau đó được gọi theo hình thức:

min_array, argmin_array = Min_Argmin(two_dim_array, axis = 0 or 1) 

Thời gian Ví dụ:

random_array = np.random.rand(20000, 20000).astype(np.float32) 

%timeit min_array, argmin_array = Min_Argmin(random_array, axis = 0) 
%timeit min_array, argmin_array = Min_Argmin(random_array, axis = 1) 

1 loops, best of 3: 405 ms per loop 
1 loops, best of 3: 307 ms per loop 

Để so sánh:

%%timeit 
min_array = random_array.min(axis = 0) 
argmin_array = random_array.argmin(axis = 0) 

1 loops, best of 3: 10.3 s per loop 

%%timeit 
min_array = random_array.min(axis = 1) 
argmin_array = random_array.argmin(axis = 1) 

1 loops, best of 3: 423 ms per loop 

Vì vậy, có một sự tăng tốc đáng kể cho axis = 0 (và vẫn còn là một lợi thế nhỏ cho axis = 1, nếu ai quan tâm đến cả tối thiểu và vị trí).