Почему einsum и NumPy быстрее, чем включает в себя встроенный в функции?
давайте начнем с трех массивов dtype=np.double. Тайминги выполняются на процессоре intel с использованием numpy 1.7.1, скомпилированного с icc и связан с intel mkl. Процессор AMD с numpy 1.6.1 скомпилирован с gcc без mkl также использовался для проверки таймингов. Обратите внимание, что тайминги масштабируются почти линейно с размером системы и не связаны с небольшими накладными расходами, понесенными в функциях numpy if заявления эти различия будут отображаться в микросекундах не миллисекунды:
arr_1D=np.arange(500,dtype=np.double)
large_arr_1D=np.arange(100000,dtype=np.double)
arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500)
arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500)
сначала давайте посмотрим на np.sum функция:
np.all(np.sum(arr_3D)==np.einsum('ijk->',arr_3D))
True
%timeit np.sum(arr_3D)
10 loops, best of 3: 142 ms per loop
%timeit np.einsum('ijk->', arr_3D)
10 loops, best of 3: 70.2 ms per loop
силы:
np.allclose(arr_3D*arr_3D*arr_3D,np.einsum('ijk,ijk,ijk->ijk',arr_3D,arr_3D,arr_3D))
True
%timeit arr_3D*arr_3D*arr_3D
1 loops, best of 3: 1.32 s per loop
%timeit np.einsum('ijk,ijk,ijk->ijk', arr_3D, arr_3D, arr_3D)
1 loops, best of 3: 694 ms per loop
внешняя сторона изделия:
np.all(np.outer(arr_1D,arr_1D)==np.einsum('i,k->ik',arr_1D,arr_1D))
True
%timeit np.outer(arr_1D, arr_1D)
1000 loops, best of 3: 411 us per loop
%timeit np.einsum('i,k->ik', arr_1D, arr_1D)
1000 loops, best of 3: 245 us per loop
все вышеперечисленное в два раза быстрее с np.einsum. Это должны быть сравнения яблок с яблоками, так как все конкретно dtype=np.double. Я ожидал бы ускорения в такой операции:
np.allclose(np.sum(arr_2D*arr_3D),np.einsum('ij,oij->',arr_2D,arr_3D))
True
%timeit np.sum(arr_2D*arr_3D)
1 loops, best of 3: 813 ms per loop
%timeit np.einsum('ij,oij->', arr_2D, arr_3D)
10 loops, best of 3: 85.1 ms per loop
Эйнсум, кажется, по крайней мере в два раза быстрее для np.inner,np.outer,np.kron, и np.sum независимо от axes выбор. Основным исключением является np.dot как он вызывает DGEMM из библиотеки BLAS. Так почему же np.einsum быстрее, чем другие функции numpy, которые эквивалентны?
случай DGEMM для полноты:
np.allclose(np.dot(arr_2D,arr_2D),np.einsum('ij,jk',arr_2D,arr_2D))
True
%timeit np.einsum('ij,jk',arr_2D,arr_2D)
10 loops, best of 3: 56.1 ms per loop
%timeit np.dot(arr_2D,arr_2D)
100 loops, best of 3: 5.17 ms per loop
ведущая теория из комментария @sebergs, что np.einsum можно использовать SSE2, но ufuncs numpy не будет до numpy 1.8 (см. изменение входа). Я считаю, что это правильный ответ, но не удалось это подтвердить. Некоторые ограниченные доказательства можно найти, изменив dtype входного массива и наблюдая разницу в скорости и тот факт, что не все наблюдают одни и те же тенденции в таймингах.
3 ответа:
во-первых, было много прошлых дискуссий об этом в списке numpy. Например, см.: http://numpy-discussion.10968.n7.nabble.com/poor-performance-of-sum-with-sub-machine-word-integer-types-td41.html http://numpy-discussion.10968.n7.nabble.com/odd-performance-of-sum-td3332.html
некоторые из них сводятся к тому, что
einsumявляется новым и, по-видимому, пытается быть лучше о выравнивании кэша и другого доступа к памяти проблемы, в то время как многие из старых функций numpy сосредоточены на легко переносимой реализации над сильно оптимизированной. Впрочем, я просто размышляю.
однако, некоторые из того, что вы делаете, не совсем сравнение "яблоки к яблокам".
в дополнение к тому, что @Jamie уже сказал,
sumиспользует более подходящий аккумулятор для массивовнапример,
sumболее тщательно о проверке типа входного сигнала и использовании соответствующий аккумулятор. Например, рассмотрим следующее:In [1]: x = 255 * np.ones(100, dtype=np.uint8) In [2]: x Out[2]: array([255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255], dtype=uint8)отметим, что
sumправильно:In [3]: x.sum() Out[3]: 25500пока
einsumдаст неверный результат:In [4]: np.einsum('i->', x) Out[4]: 156но если мы используем менее ограничен
dtype, мы все равно получим результат, который вы ожидали бы:In [5]: y = 255 * np.ones(100) In [6]: np.einsum('i->', y) Out[6]: 25500.0
теперь, когда numpy 1.8 выпущен, где согласно документам все ufuncs должны использовать SSE2, я хотел дважды проверить, что комментарий Себерга о SSE2 был действителен.
для выполнения теста была создана новая установка python 2.7 - numpy 1.7 и 1.8 были скомпилированы с
iccиспользование стандартных опций на ядре AMD opteron под управлением Ubuntu.это тестовый запуск как до, так и после обновления 1.8:
import numpy as np import timeit arr_1D=np.arange(5000,dtype=np.double) arr_2D=np.arange(500**2,dtype=np.double).reshape(500,500) arr_3D=np.arange(500**3,dtype=np.double).reshape(500,500,500) print 'Summation test:' print timeit.timeit('np.sum(arr_3D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print timeit.timeit('np.einsum("ijk->", arr_3D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print '----------------------\n' print 'Power test:' print timeit.timeit('arr_3D*arr_3D*arr_3D', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print timeit.timeit('np.einsum("ijk,ijk,ijk->ijk", arr_3D, arr_3D, arr_3D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print '----------------------\n' print 'Outer test:' print timeit.timeit('np.outer(arr_1D, arr_1D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print timeit.timeit('np.einsum("i,k->ik", arr_1D, arr_1D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print '----------------------\n' print 'Einsum test:' print timeit.timeit('np.sum(arr_2D*arr_3D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print timeit.timeit('np.einsum("ij,oij->", arr_2D, arr_3D)', 'import numpy as np; from __main__ import arr_1D, arr_2D, arr_3D', number=5)/5 print '----------------------\n'Numpy 1.7.1:
Summation test: 0.172988510132 0.0934836149216 ---------------------- Power test: 1.93524689674 0.839519000053 ---------------------- Outer test: 0.130380821228 0.121401786804 ---------------------- Einsum test: 0.979052495956 0.126066613197Numpy 1.8:
Summation test: 0.116551589966 0.0920487880707 ---------------------- Power test: 1.23683619499 0.815982818604 ---------------------- Outer test: 0.131808176041 0.127472200394 ---------------------- Einsum test: 0.781750011444 0.129271841049Я думаю, что это довольно убедительно, что SSE играет большую роль в разнице во времени, следует отметить, что повторение этих тестов тайминги очень только ~0.003 С. оставшаяся разница должна быть покрыта в других ответах на этот вопрос.
Я думаю, что эти тайминги объяснить, что происходит:
a = np.arange(1000, dtype=np.double) %timeit np.einsum('i->', a) 100000 loops, best of 3: 3.32 us per loop %timeit np.sum(a) 100000 loops, best of 3: 6.84 us per loop a = np.arange(10000, dtype=np.double) %timeit np.einsum('i->', a) 100000 loops, best of 3: 12.6 us per loop %timeit np.sum(a) 100000 loops, best of 3: 16.5 us per loop a = np.arange(100000, dtype=np.double) %timeit np.einsum('i->', a) 10000 loops, best of 3: 103 us per loop %timeit np.sum(a) 10000 loops, best of 3: 109 us per loopтаким образом, у вас в основном есть почти постоянные накладные расходы 3us при вызове
np.sumovernp.einsum, так что они в основном работают так же быстро, но один занимает немного больше времени, чтобы идти. Почему это может быть? Мои деньги на следующее:a = np.arange(1000, dtype=object) %timeit np.einsum('i->', a) Traceback (most recent call last): ... TypeError: invalid data type for einsum %timeit np.sum(a) 10000 loops, best of 3: 20.3 us per loopне уверен, что именно происходит, но кажется, что
np.einsumпропускает некоторые проверки для извлечения конкретных функций типа для выполнения умножения и сложения, а также идем прямо с*и+для стандартных типов C только.
многомерные случаи не отличаются:
n = 10; a = np.arange(n**3, dtype=np.double).reshape(n, n, n) %timeit np.einsum('ijk->', a) 100000 loops, best of 3: 3.79 us per loop %timeit np.sum(a) 100000 loops, best of 3: 7.33 us per loop n = 100; a = np.arange(n**3, dtype=np.double).reshape(n, n, n) %timeit np.einsum('ijk->', a) 1000 loops, best of 3: 1.2 ms per loop %timeit np.sum(a) 1000 loops, best of 3: 1.23 ms per loopтаким образом, в основном постоянные накладные расходы, а не более быстрый ход, как только они доберутся до него.