Ответ 1
Сначала я подозревал, что numpy вызывает BLAS, но, по крайней мере, на моей машине (python 2.7.13, numpy 1.11.2, OpenBLAS), это не так, поскольку быстрая проверка с помощью gdb показывает:
> gdb --args python timing.py
...
Size: 100000000.0
^C
Thread 1 "python" received signal SIGINT, Interrupt.
sse2_binary_scalar2_divide_DOUBLE (op=0x7fffb3aee010, ip1=0x7fffb3aee010, ip2=0x6fe2c0, n=100000000)
at numpy/core/src/umath/simd.inc.src:491
491 numpy/core/src/umath/simd.inc.src: No such file or directory.
(gdb) disass
...
0x00007fffe6ea6228 <+392>: movapd (%rsi,%rax,8),%xmm0
0x00007fffe6ea622d <+397>: divpd %xmm1,%xmm0
=> 0x00007fffe6ea6231 <+401>: movapd %xmm0,(%rdi,%rax,8)
...
(gdb) p $xmm1
$1 = {..., v2_double = {0.5, 0.5}, ...}
Фактически, numpy выполняет точно такой же общий цикл независимо от используемой константы. Таким образом, все временные разности обусловлены исключительно процессором.
Фактически, деление - это команда с очень переменным временем выполнения. Объем выполняемой работы зависит от битовых паттернов операндов, а также могут быть обнаружены и ускорены специальные случаи. Согласно эти таблицы (точность которых я не знаю), на вашем E5-2620 (Sandy Bridge) DIVPD имеет задержку и обратную пропускную способность из 10-22 циклов, а MULPS имеет латентность 10 циклов и обратную пропускную способность 5 циклов.
Теперь, если A*2.0 медленнее, чем A*=2.0. gdb показывает, что для умножения используется точно такая же функция, за исключением того, что вывод op отличается от первого входа ip1. Таким образом, это должно быть просто артефактом дополнительной памяти, которая втягивается в кеш, замедляя операцию без входа для большого ввода (хотя MULPS производит только 2 * 8/5 = 3,2 байта вывода за цикл!). При использовании буферов размером 1е4 все вписывается в кеш, что не оказывает существенного эффекта, а другие накладные расходы в основном заглушают разницу между A/=0.5 и A/=0.51.
Тем не менее, в этих таймингах есть много странных эффектов, поэтому я построил график (код для его создания ниже)
Я построил размер массива A против количества циклов ЦП в инструкции DIVPD/MULPD/ADDPD. Я побежал на 3,3 ГГц AMD FX-6100. Желтые и красные вертикальные линии - это размер кеша L2 и L3. Синяя линия - это предполагаемая максимальная пропускная способность DIVPD в соответствии с этими таблицами, 1/4.5cycles (что кажется сомнительным). Как вы можете видеть, даже не A+=2.0 приближается к этому, даже когда "Накладные расходы" выполнения операции numpy приближаются к нулю. Таким образом, существует около 24 циклов накладных расходов, только цикл и чтение и запись 16 байт в/из кеша L2! Довольно шокирующий, возможно, доступ к памяти не согласован.
Множество интересных эффектов:
- Ниже массивов 30 Кбайт большая часть времени является накладными в python/numpy
- Умножение и добавление имеют одинаковую скорость (как указано в таблицах Agner).
- Разница в скорости между
A/=0.5иA/=0.51падает вправо от графика; это связано с тем, что, когда время чтения/записи увеличивается, оно перекрывается и маскирует некоторое время, затрачиваемое на деление. По этой причинеA/=0.5,A*=2.0иA+=2.0становятся равной скорости. - Сравнивая максимальную разницу между
A/=0.51,A/=0.5иA+=2.0, предполагает, что деление имеет пропускную способность 4,5-44 цикла, что не соответствует 4,5-11 в таблице Agner. - Однако разница между A/= 0,5 и A/= 0,51 в основном исчезает, когда накладные расходы numpy становятся большими, хотя разница в нескольких циклах еще меньше. Это трудно объяснить, потому что избыточные служебные данные не могут маскировать время для деления.
- Операции, которые не являются на месте (пунктирные линии), становятся невероятно медленными, когда они намного больше, чем размер кеша L3, но операции на месте этого не делают. Они требуют удвоить пропускную способность памяти в ОЗУ, но я не могу объяснить, почему они будут в 20 раз медленнее!
- Пунктирные линии расходятся слева. Это предположительно, потому что деление и умножение обрабатываются различными функциями numpy с разным объемом служебных данных.
К сожалению, на другой машине с процессором с разной скоростью FPU, размером кеша, пропускной способностью памяти, версией numpy и т.д. эти кривые могут выглядеть совсем по-другому.
Мой взлет от этого: объединение нескольких арифметических операций с numpy будет во много раз медленнее, чем делать то же самое в Cython, итерации по входам только один раз, потому что нет "сладкого пятна", при котором стоимость арифметических операций доминирует над другими затратами.
import numpy as np
import timeit
import matplotlib.pyplot as plt
CPUHz = 3.3e9
divpd_cycles = 4.5
L2cachesize = 2*2**20
L3cachesize = 8*2**20
def timeit_command(command, pieces, size):
return min(timeit.repeat("for i in xrange(%d): %s" % (pieces, command),
"import numpy; A = numpy.random.rand(%d)" % size, number = 6))
def run():
totaliterations = 1e7
commands=["A/=0.5", "A/=0.51", "A/0.5", "A*=2.0", "A*2.0", "A+=2.0"]
styles=['-', '-', '--', '-', '--', '-']
def draw_graph(command, style, compute_overhead = False):
sizes = []
y = []
for pieces in np.logspace(0, 5, 11):
size = int(totaliterations / pieces)
sizes.append(size * 8) # 8 bytes per double
time = timeit_command(command, pieces, (4 if compute_overhead else size))
# Divide by 2 because SSE instructions process two doubles each
cycles = time * CPUHz / (size * pieces / 2)
y.append(cycles)
if compute_overhead:
command = "numpy overhead"
plt.semilogx(sizes, y, style, label = command, linewidth = 2, basex = 10)
plt.figure()
for command, style in zip(commands, styles):
print command
draw_graph(command, style)
# Plot overhead
draw_graph("A+=1.0", '-', compute_overhead=True)
plt.legend(loc = 'best', prop = {'size':9}, handlelength = 3)
plt.xlabel('Array size in bytes')
plt.ylabel('CPU cycles per SSE instruction')
# Draw vertical and horizontal lines
ymin, ymax = plt.ylim()
plt.vlines(L2cachesize, ymin, ymax, color = 'orange', linewidth = 2)
plt.vlines(L3cachesize, ymin, ymax, color = 'red', linewidth = 2)
xmin, xmax = plt.xlim()
plt.hlines(divpd_cycles, xmin, xmax, color = 'blue', linewidth = 2)