Ваше предположение, что ускорение происходит только / главным образом из-за распараллеливания, неверно.Как уже указывал @Brenlla, львиная доля ускорения Numberxpr обычно достигается за счет лучшего использования кэша.Однако есть и другие причины.
Во-первых, numpy и Numberxpr не оценивают одно и то же выражение:
- numpy вычисляет
x**3 и x**2 как pow(x,3)и pow(x,2). - Numberxpr позволяет себе оценить его как
x**3=x*x*x и x**2=x*x.
pow сложнее, чем одно или два умножения и, следовательно, намногомедленнее, сравните:
ne.set_num_threads(1)
%timeit ne.evaluate("0.25*x**3 + 0.75*x**2 + 1.5*x - 2")
# 60.7 ms ± 1.2 ms, base line on my machine
%timeit 0.25*x**3 + 0.75*x**2 + 1.5*x - 2
# 766 ms ± 4.02 ms
%timeit 0.25*x*x*x + 0.75*x*x + 1.5*x - 2
# 130 ms ± 692 µs
Теперь, Numberxpr только в два раза быстрее.Я предполагаю, что версия pow была привязана к процессору, а версия умножения больше связана с памятью.
Numexpr в основном светит, когда данные большие - больше, чем кэш L3 (например, 15 МБ на моем компьютере), что указано в вашем примере, так как x составляет около 76 МБ:
- Numberxp оценивает по блокам - т.е. все операции оцениваются для блока, и каждый блок подходит(по крайней мере) L3-кеш, таким образом максимизируя использование кеша.ТОЛЬКО после завершения одного блока выполняется оценка другого блока.
- numpy оценивает одну операцию за другой для всех данных, таким образом, данные удаляются из кэша, прежде чем их можно будет повторно использовать.
Мы можем посмотреть на ошибки в кеше, используя, например, valgrind (см. Сценарии в приложении к этому сообщению):
>>> valgrind --tool=cachegrind python np_version.py
...
...
==5676== D refs: 1,144,572,370 (754,717,376 rd + 389,854,994 wr)
==5676== D1 misses: 220,844,716 (181,436,970 rd + 39,407,746 wr)
==5676== LLd misses: 217,056,340 (178,062,890 rd + 38,993,450 wr)
==5676== D1 miss rate: 19.3% ( 24.0% + 10.1% )
==5676== LLd miss rate: 19.0% ( 23.6% + 10.0% )
....
Интересная часть для нас - LLd-misses (то есть ошибки в L3,см. здесь для получения информации о интерпретации выходных данных) - около 25% обращений к чтению являются пропущенными.
Тот же анализ для numbersxpr показывает:
>>> valgrind --tool=cachegrind python ne_version.py
...
==5145== D refs: 2,612,495,487 (1,737,673,018 rd + 874,822,469 wr)
==5145== D1 misses: 110,971,378 ( 86,949,951 rd + 24,021,427 wr)
==5145== LLd misses: 29,574,847 ( 15,579,163 rd + 13,995,684 wr)
==5145== D1 miss rate: 4.2% ( 5.0% + 2.7% )
==5145== LLd miss rate: 1.1% ( 0.9% + 1.6% )
...
Только 5% операций чтения - это промахи!
Однако у numpy есть некоторые преимущества: под капотом numpy использует mkl-рутины (по крайней мере, на моей машине), а Numberxpr - нет.Таким образом, numpy заканчивается использованием упакованных SSE-операций (movups + mulpd + addpd), в то время как Numberxpr заканчивается скалярными версиями (movsd + mulsd).
Это объясняет 25-процентную вероятность пропусков в numpy-версии: одно чтение составляет 128 бит (movups), что означает, что после 4 операций чтения строка кэша (64 байта) обрабатывается и создается ошибка.Это можно увидеть в профиле (например, perf в Linux):
32,93 │ movups 0x10(%r15,%rcx,8),%xmm4
1,33 │ movups 0x20(%r15,%rcx,8),%xmm5
1,71 │ movups 0x30(%r15,%rcx,8),%xmm6
0,76 │ movups 0x40(%r15,%rcx,8),%xmm7
24,68 │ movups 0x50(%r15,%rcx,8),%xmm8
1,21 │ movups 0x60(%r15,%rcx,8),%xmm9
2,54 │ movups 0x70(%r15,%rcx,8),%xmm10
каждый четвертый movups требует больше времени, потому что он ожидает доступа к памяти.
Numpy сияет для меньших размеров массива, которые соответствуют кэш-памяти L1 (но тогда львиная доля приходится на накладные расходы, а не на сами вычисления, которые быстрее в numpy - но это не играет большой роли):
x = np.linspace(-1, 1, 10**3)
%timeit ne.evaluate("0.25*x*x*x + 0.75*x*x + 1.5*x - 2")
# 20.1 µs ± 306 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit 0.25*x*x*x + 0.75*x*x + 1.5*x - 2
# 13.1 µs ± 125 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
Примечание: было бы быстрее оценить функцию как ((0.25*x + 0.75)*x + 1.5)*x - 2.
Оба из-за меньшей загрузки ЦП:
# small x - CPU bound
x = np.linspace(-1, 1, 10**3)
%timeit ((0.25*x + 0.75)*x + 1.5)*x - 2
# 9.02 µs ± 204 ns
и меньше обращений к памяти:
# large x - memory bound
x = np.linspace(-1, 1, 10**7)
%timeit ((0.25*x + 0.75)*x + 1.5)*x - 2
# 73.8 ms ± 3.71 ms
Списки:
A np_version.py:
import numpy as np
x = np.linspace(-1, 1, 10**7)
for _ in range(10):
cubic_poly = 0.25*x*x*x + 0.75*x*x + 1.5*x - 2
B ne_version.py:
import numpy as np
import numexpr as ne
x = np.linspace(-1, 1, 10**7)
ne.set_num_threads(1)
for _ in range(10):
ne.evaluate("0.25*x**3 + 0.75*x**2 + 1.5*x - 2")