Почему numpy Крон так быстро? - PullRequest
Новая
       57

Почему numpy Крон так быстро?

1 голос
/ 24 февраля 2020

Я пытался реализовать функцию kronecker product . Ниже приведены три идеи, которые у меня есть: 

def kron(arr1, arr2):
    """columnwise outer product, avoiding relocate elements.
    """
    r1, c1 = arr1.shape
    r2, c2 = arr2.shape
    nrows, ncols = r1 * r2, c1 * c2
    res = np.empty((nrows, ncols))
    for idx1 in range(c1):
        for idx2 in range(c2):
            new_c = idx1 * c2 + idx2
            temp = np.zeros((r2, r1))
            temp_kron = scipy.linalg.blas.dger(
                alpha=1.0, x=arr2[:, idx2], y=arr1[:, idx1], incx=1, incy=1, 
                a=temp)
            res[:, new_c] = np.ravel(temp_kron, order='F')
    return res

def kron2(arr1, arr2):
    """First outer product, then rearrange items.
    """
    r1, c1 = arr1.shape
    r2, c2 = arr2.shape
    nrows, ncols = r1 * r2, c1 * c2
    tmp = np.outer(arr2, arr1)
    res = np.empty((nrows, ncols))
    for idx in range(arr1.size):
        for offset in range(c2):
            orig = tmp[offset::c2, idx]
            dest_coffset = idx % c1 * c2 + offset
            dest_roffset = (idx // c1) * r2
            res[dest_roffset:dest_roffset+r2, dest_coffset] = orig
    return res

def kron3(arr1, arr2):
    """First outer product, then rearrange items.
    """
    r1, c1 = arr1.shape
    r2, c2 = arr2.shape
    nrows, ncols = r1 * r2, c1 * c2
    tmp = np.outer(np.ravel(arr2, 'F'), np.ravel(arr1, 'F'))
    res = np.empty((nrows, ncols))
    for idx in range(arr1.size):
        for offset in range(c2):
            orig_offset = offset * r2
            orig = tmp[orig_offset:orig_offset+r2, idx]
            dest_c = idx // r1 * c2 + offset
            dest_r = idx % r1 * r2
            res[dest_r:dest_r+r2, dest_c] = orig
    return res

На основе этой записи stackoverflow я создал декоратор MeasureTime. Естественным эталоном будет сравнение с numpy .kron . Ниже приведены мои тестовые функции: 

@MeasureTime
def test_np_kron(arr1, arr2, number=1000):
    for _ in range(number):
        np.kron(arr1, arr2)
    return

@MeasureTime
def test_kron(arr1, arr2, number=1000):
    for _ in range(number):
        kron(arr1, arr2)

@MeasureTime
def test_kron2(arr1, arr2, number=1000):
    for _ in range(number):
        kron2(arr2, arr1)

@MeasureTime
def test_kron3(arr1, arr2, number=1000):
    for _ in range(number):
        kron3(arr2, arr1)

Оказалось, что функции крон Numpy намного лучше: 

arr1 = np.array([[1,-4,7], [-2, 3, 3]], dtype=np.float64, order='F')
arr2 = np.array([[8, -9, -6, 5], [1, -3, -4, 7], [2, 8, -8, -3], [1, 2, -5, -1]], dtype=np.float64, order='F')

In [243]: test_np_kron(arr1, arr2, number=10000)
Out [243]: "test_np_kron": 0.19688990000577178s

In [244]: test_kron(arr1, arr2, number=10000)
Out [244]: "test_kron": 0.6094115000014426s

In [245]: test_kron2(arr1, arr2, number=10000)
Out [245]: "test_kron2": 0.5699560000066413s

In [246]: test_kron3(arr1, arr2, number=10000)
Out [246]: "test_kron3": 0.7134822000080021s

Я хотел бы знать, почему это так? Это потому, что метод изменения формы Numpy гораздо эффективнее, чем ручное копирование (хотя по-прежнему используется numpy)? Я был озадачен, так как в противном случае я также использовал np.outer / blas.dger. Единственное отличие, которое я здесь заметил, заключалось в том, как мы организовали конечные результаты. Почему перефразировка NumPy делает это хорошо?

Здесь - это ссылка на NumPy 1,17 крон источника.

Обновления: Забыл упомянуть, во-первых, что я пытался создать прототип в python, а затем внедрить kron с использованием C ++ с помощью cblas / lapack. Имел ли какой-то существующий «крон» необходимость в рефакторинге. Затем я наткнулся на Numpy изменить форму и был действительно впечатлен.

Заранее спасибо за ваше время!

1 Ответ

3 голосов
/ 25 февраля 2020

Давайте поэкспериментируем с 2 маленькими массивами: 

In [124]: A, B = np.array([[1,2],[3,4]]), np.array([[10,11],[12,13]])

kron производит: 

In [125]: np.kron(A,B)                                                                         
Out[125]: 
array([[10, 11, 20, 22],
       [12, 13, 24, 26],
       [30, 33, 40, 44],
       [36, 39, 48, 52]])

outer производит те же числа, но с другим расположением: 

In [126]: np.outer(A,B)                                                                        
Out[126]: 
array([[10, 11, 12, 13],
       [20, 22, 24, 26],
       [30, 33, 36, 39],
       [40, 44, 48, 52]])

kron преобразует его в комбинацию форм A и B: 

In [127]: np.outer(A,B).reshape(2,2,2,2)                                                       
Out[127]: 
array([[[[10, 11],
         [12, 13]],

        [[20, 22],
         [24, 26]]],


       [[[30, 33],
         [36, 39]],

        [[40, 44],
         [48, 52]]]])

, затем он объединяет 4 измерения в 2 с concatenate: 

In [128]: np.concatenate(np.concatenate(_127, 1),1)                                            
Out[128]: 
array([[10, 11, 20, 22],
       [12, 13, 24, 26],
       [30, 33, 40, 44],
       [36, 39, 48, 52]])

Альтернативой является замена осей и изменение формы: 

In [129]: _127.transpose(0,2,1,3).reshape(4,4)                                                 
Out[129]: 
array([[10, 11, 20, 22],
       [12, 13, 24, 26],
       [30, 33, 40, 44],
       [36, 39, 48, 52]])

При первом изменении формы и транспонировании создается представление, но при втором изменении формы создается копия. Конкатенация делает копию. Но все эти действия выполняются в скомпилированном numpy коде.

Определение функций: 

def foo1(A,B):
    temp = np.outer(A,B)
    temp = temp.reshape(A.shape + B.shape)
    return np.concatenate(np.concatenate(temp, 1), 1)

def foo2(A,B):
    temp = np.outer(A,B)
    nz = temp.shape
    temp = temp.reshape(A.shape + B.shape)
    return temp.transpose(0,2,1,3).reshape(nz)

тестирование: 

In [141]: np.allclose(np.kron(A,B), foo1(A,B))                                                 
Out[141]: True
In [142]: np.allclose(np.kron(A,B), foo2(A,B))                                                 
Out[142]: True

синхронизация: 

In [143]: timeit np.kron(A,B)                                                                  
42.4 µs ± 294 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [145]: timeit foo1(A,B)                                                                     
26.3 µs ± 38.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [146]: timeit foo2(A,B)                                                                     
13.8 µs ± 19.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Мой код может нуждаться в некотором обобщении, но это демонстрирует правильность подхода.

===

С вашим kron: 

In [150]: kron(A,B)                                                                            
Out[150]: 
array([[10., 11., 20., 22.],
       [12., 13., 24., 26.],
       [30., 33., 40., 44.],
       [36., 39., 48., 52.]])
In [151]: timeit kron(A,B)                                                                     
55.3 µs ± 1.59 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
...