Почему функция numpy einsum медленнее, чем встроенные функции numpy?

Обычно я получал хорошую производительность из функции einsum numpy (и мне нравится ее синтаксис). @ Ответ Офиона на этот вопрос показывает, что – для проверенных случаев – einsum последовательно превосходит «встроенные» функции (иногда немного, иногда по большому счету). Но я просто столкнулся с ситуацией, когда einsum намного медленнее. Рассмотрим следующие эквивалентные функции:

(M, K) = (1000000, 20) C = np.random.rand(K, K) X = np.random.rand(M, K) def func_dot(C, X): Y = X.dot(C) return np.sum(Y * X, axis=1) def func_einsum(C, X): return np.einsum('ik,km,im->i', X, C, X) def func_einsum2(C, X): # Like func_einsum but break it into two steps. A = np.einsum('ik,km', X, C) return np.einsum('ik,ik->i', A, X) 

Я ожидал, что func_einsum будет работать быстрее, но это не то, с чем я сталкиваюсь. Запуск на четырехъядерном процессоре с гиперпотоком, numpy версии 1.9.0.dev-7ae0206 и многопоточность с OpenBLAS, я получаю следующие результаты:

  • «Клонирование» строк или столбцов векторов
  • Как найти линейно независимые строки из матрицы
  • Эффективная и pythonic проверка для сингулярной матрицы
  • Python Обратная матрица
  • Вычислить нормальную форму матрицы Джорда в Python / NumPy
  • Векторизация этого двойного цикла Numpy
  •  In [2]: %time y1 = func_dot(C, X) CPU times: user 320 ms, sys: 312 ms, total: 632 ms Wall time: 209 ms In [3]: %time y2 = func_einsum(C, X) CPU times: user 844 ms, sys: 0 ns, total: 844 ms Wall time: 842 ms In [4]: %time y3 = func_einsum2(C, X) CPU times: user 292 ms, sys: 44 ms, total: 336 ms Wall time: 334 ms 

    Когда я увеличиваю K до 200, различия более экстремальные:

     In [2]: %time y1= func_dot(C, X) CPU times: user 4.5 s, sys: 1.02 s, total: 5.52 s Wall time: 2.3 s In [3]: %time y2= func_einsum(C, X) CPU times: user 1min 16s, sys: 44 ms, total: 1min 16s Wall time: 1min 16s In [4]: %time y3 = func_einsum2(C, X) CPU times: user 15.3 s, sys: 312 ms, total: 15.6 s Wall time: 15.6 s 

    Может ли кто-нибудь объяснить, почему einsum здесь намного медленнее?

    Если это имеет значение, вот моя конфигурация numpy:

     In [6]: np.show_config() lapack_info: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] language = f77 atlas_threads_info: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('ATLAS_WITHOUT_LAPACK', None)] language = c include_dirs = ['/usr/local/include'] blas_opt_info: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('ATLAS_INFO', '"\\"None\\""')] language = c include_dirs = ['/usr/local/include'] atlas_blas_threads_info: libraries = ['openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('ATLAS_INFO', '"\\"None\\""')] language = c include_dirs = ['/usr/local/include'] lapack_opt_info: libraries = ['openblas', 'openblas'] library_dirs = ['/usr/local/lib'] define_macros = [('ATLAS_WITHOUT_LAPACK', None)] language = f77 include_dirs = ['/usr/local/include'] lapack_mkl_info: NOT AVAILABLE blas_mkl_info: NOT AVAILABLE mkl_info: NOT AVAILABLE 

  • Поиск индекса ближайшей точки в массивах numpy координат x и y
  • Как вернуть все минимальные индексы в numpy
  • Scipy интерполяция, как изменить размер / перемасштабировать матрицу 3x3 до 5x5?
  • Невозможно использовать 128-битный float в Python на 64-битной архитектуре
  • Matplotlib: ValueError: x и y должны иметь одинаковое первое измерение
  • StringIO в python3
  • 2 Solutions collect form web for “Почему функция numpy einsum медленнее, чем встроенные функции numpy?”

    Вы можете получить лучшее из обоих миров:

     def func_dot_einsum(C, X): Y = X.dot(C) return np.einsum('ij,ij->i', Y, X) 

    В моей системе:

     In [7]: %timeit func_dot(C, X) 10 loops, best of 3: 31.1 ms per loop In [8]: %timeit func_einsum(C, X) 10 loops, best of 3: 105 ms per loop In [9]: %timeit func_einsum2(C, X) 10 loops, best of 3: 43.5 ms per loop In [10]: %timeit func_dot_einsum(C, X) 10 loops, best of 3: 21 ms per loop 

    Когда доступно, np.dot использует BLAS, MKL или любую np.dot вас библиотеку. Таким образом, вызов np.dot почти наверняка многопоточен. np.einsum имеет свои собственные циклы, поэтому не использует ни одну из этих оптимизаций, кроме собственного использования SIMD, чтобы ускорить работу над реализацией ванили C.


    Затем есть многопользовательский вызов einsum, который работает намного медленнее … Источник numpy для einsum очень сложный, и я не совсем понимаю его. Поэтому имейте в виду, что в лучшем случае спекулятивно, но вот что я думаю, происходит …

    Когда вы запускаете что-то вроде np.einsum('ij,ij->i', a, b) , преимущество над выполнением np.sum(a*b, axis=1) происходит из-за того, что не нужно создавать промежуточный массив со всеми продуктов и зацикливания над ним. Так что на низком уровне все происходит так:

     for i in range(I): out[i] = 0 for j in range(J): out[i] += a[i, j] * b[i, j] 

    Скажите, что вы после чего-то вроде:

     np.einsum('ij,jk,ik->i', a, b, c) 

    Вы можете выполнить ту же операцию, что и

     np.sum(a[:, :, None] * b[None, :, :] * c[:, None, :], axis=(1, 2)) 

    И я думаю, что einsum делает это, чтобы запустить этот последний код без необходимости создания экземпляра огромного промежуточного массива, что, безусловно, ускоряет работу:

     In [29]: a, b, c = np.random.rand(3, 100, 100) In [30]: %timeit np.einsum('ij,jk,ik->i', a, b, c) 100 loops, best of 3: 2.41 ms per loop In [31]: %timeit np.sum(a[:, :, None] * b[None, :, :] * c[:, None, :], axis=(1, 2)) 100 loops, best of 3: 12.3 ms per loop 

    Но если вы внимательно посмотрите на это, избавление от промежуточного хранения может быть ужасным. Это то, что я думаю, что einsum делает на низком уровне:

     for i in range(I): out[i] = 0 for j in range(J): for k in range(K): out[i] += a[i, j] * b[j, k] * c[i, k] 

    Но вы повторяете массу операций! Если вы это сделали:

     for i in range(I): out[i] = 0 for j in range(J): temp = 0 for k in range(K): temp += b[j, k] * c[i, k] out[i] += a[i, j] * temp 

    вы бы делали I * J * (K-1) меньше умножений (и дополнительных дополнений I * J ), и сэкономьте себе массу времени. Я предполагаю, что einsum недостаточно умен, чтобы оптимизировать ситуацию на этом уровне. В исходном коде есть намек на то, что он только оптимизирует операции с 1 или 2 операндами, а не 3. В любом случае автоматизация этого для общих входов кажется чем-то простым, но простым …

    einsum имеет специализированный случай для '2 операндов, ndim = 2'. В этом случае есть 3 операнда и всего 3 измерения. Поэтому он должен использовать общий nditer .

    При попытке понять, как разбирается ввод строки, я написал чистый симулятор Pinson einsum, https://github.com/hpaulj/numpy-einsum/blob/master/einsum_py.py

    Функции einsum и sum-of-products (урезанные):

     def myeinsum(subscripts, *ops, **kwargs): # dropin preplacement for np.einsum (more or less) <parse subscript strings> <prepare op_axes> x = sum_of_prod(ops, op_axes, **kwargs) return x def sum_of_prod(ops, op_axes,...): ... it = np.nditer(ops, flags, op_flags, op_axes) it.operands[nop][...] = 0 it.reset() for (x,y,z,w) in it: w[...] += x*y*z return it.operands[nop] 

    myeinsum('ik,km,im->i',X,C,X,debug=True) вывод для myeinsum('ik,km,im->i',X,C,X,debug=True) с (M,K)=(10,5)

     {'max_label': 109, 'min_label': 105, 'nop': 3, 'shapes': [(10, 5), (5, 5), (10, 5)], ....}} ... iter labels: [105, 107, 109],'ikm' op_axes [[0, 1, -1], [-1, 0, 1], [0, -1, 1], [0, -1, -1]] 

    Если вы напишете функцию sum-of-prod подобную этой в cython вы должны получить что-то близкое к обобщенному einsum .

    С полным (M,K) этот симулированный einsum на 6-7 раз медленнее.


    Некоторые тайминги, основанные на других ответах:

     In [84]: timeit np.dot(X,C) 1 loops, best of 3: 781 ms per loop In [85]: timeit np.einsum('ik,km->im',X,C) 1 loops, best of 3: 1.28 s per loop In [86]: timeit np.einsum('im,im->i',A,X) 10 loops, best of 3: 163 ms per loop 

    Этот 'im,im->i' step is substantially faster than the other. The sum dimension, 'im,im->i' step is substantially faster than the other. The sum dimension, m is only 20. I suspect einsum` рассматривает это как частный случай.

     In [87]: timeit np.einsum('im,im->i',np.dot(X,C),X) 1 loops, best of 3: 950 ms per loop In [88]: timeit np.einsum('im,im->i',np.einsum('ik,km->im',X,C),X) 1 loops, best of 3: 1.45 s per loop 

    Время для этих составных вычислений представляет собой просто суммы соответствующих кусков.

    Python - лучший язык программирования в мире.