numpy - einsum vs наивная реализация реализована - PullRequest
Новая
       5

numpy - einsum vs наивная реализация реализована

0 голосов
/ 18 марта 2020

У меня есть двумерный массив Y размера (N, M), скажем, например: 

N, M = 200, 100
Y = np.random.normal(0,1,(N,M))

Для каждого N я хочу вычислить скалярное произведение вектора (M, 1) с его транспонированием, которое возвращает (M, M) матрицу. Один из способов сделать это неэффективно: 

Y = Y[:,:,np.newaxis]
[Y[i,:,:] @ Y[i,:,:].T for i in range(N)]

, что довольно медленно: timeit во второй строке возвращает 

11.7 ms ± 1.39 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

Я думал, что гораздо лучший способ сделать это - использовать Функция einsum numpy (https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html): 

np.einsum('ijk,imk->ijm', Y, Y, optimize=True)

(что означает: для каждой строки i создайте (j, k) матрицу, в которой ее элементы получаются из точки произведение в последнем измерении m)

Два метода возвращают один и тот же результат, но время выполнения этой новой версии разочаровывает (только чуть более чем в два раза быстрее) 

3.82 ms ± 146 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Можно ожидать гораздо большего улучшения с использованием векторизованной функции einsum, поскольку первый метод очень неэффективен ... У вас есть объяснение этому? Существует ли лучший способ сделать этот расчет?

Ответы [ 2 ]

1 голос
/ 18 марта 2020 
In [60]: N, M = 200, 100 
    ...: Y = np.random.normal(0,1,(N,M))                                                                             
In [61]: Y1 = Y[:,:,None]

Ваша итерация, 200 шагов для создания (100 100) массивов: 

In [62]: timeit [Y1[i,:,:]@Y1[i,:,:].T for i in range(N)]                                                            
18.5 ms ± 784 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

einsum только немного быстрее: 

In [64]: timeit np.einsum('ijk,imk->ijm', Y1,Y1)                                                                     
14.5 ms ± 114 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

, но вы можете применить @ в полном «пакетном» режиме с: 

In [65]: timeit Y[:,:,None]@Y[:,None,:]                                                                              
7.63 ms ± 224 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Но, как отмечает Дивакар, ось суммы имеет размер 1, так что вы можете использовать обычное умноженное вещание. Это внешний продукт, а не матричный. 

In [66]: timeit Y[:,:,None]*Y[:,None,:]                                                                              
8.2 ms ± 64.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

«векторизация» дает большие выгоды при выполнении множества итераций простой операции. Для меньшего количества операций в более сложных операциях усиление не так велико.

0 голосов
/ 19 марта 2020

Это старый пост, но он охватывает тему во многих деталях: эффективный внешний продукт .

В частности, если вы заинтересованы в добавлении зависимости от numba, это может быть вашим самым быстрым вариантом .

Обновление части кода numba из оригинального поста и добавление мульти-внешнего продукта: 

import numpy as np
from numba import jit
from numba.typed import List

@jit(nopython=True)
def outer_numba(a, b):
    m = a.shape[0]
    n = b.shape[0]
    result = np.empty((m, n))
    for i in range(m):
        for j in range(n):
            result[i, j] = a[i]*b[j]
    return result

@jit(nopython=True)
def multi_outer_numba(Y):

    all_result = List()
    for k in range(Y.shape[0]):
        y = Y[k]           
        n = y.shape[0]
        tmp_res = np.empty((n, n))
        for i in range(n):
            for j in range(n):
                tmp_res[i, j] = y[i]*y[j]
        all_result.append(tmp_res)

    return all_result

r = [outer_numba(Y[i],Y[i]) for i in range(N)]
r = multi_outer_numba(Y)
...