Подтвердить что ты не робот

Numpy: Что особенного в делении на 0,5?

В этом ответе @Dunes говорится, что из-за конвейерности (почти) нет разницы между умножением и делением с плавающей запятой. Однако из моего опыта с другими языками я ожидал бы, что деление будет медленнее.

Мой маленький тест выглядит следующим образом:

A=np.random.rand(size)
command(A)

Для разных команд и size=1e8 я получаю следующие моменты на моей машине:

Command:    Time[in sec]:
A/=0.5      2.88435101509
A/=0.51     5.22591209412
A*=2.0      1.1831600666
A*2.0       3.44263911247  //not in-place, more cache misses?
A+=A        1.2827270031

Самая интересная часть: деление на 0.5 почти в два раза быстрее, чем деление на 0.51. Можно предположить, что это связано с некоторой интеллектуальной оптимизацией, например. заменяя деление на A+A. Однако таймеры A*2 и A+A слишком далеко, чтобы поддержать это утверждение.

В общем случае деление по поплавкам со значениями (1/2)^n выполняется быстрее:

Size: 1e8
    Command:    Time[in sec]:
    A/=0.5      2.85750007629
    A/=0.25     2.91607499123
    A/=0.125    2.89376401901
    A/=2.0      2.84901714325
    A/=4.0      2.84493684769
    A/=3.0      5.00480890274
    A/=0.75     5.0354950428
    A/=0.51     5.05687212944

Это становится еще интереснее, если мы посмотрим на size=1e4:

Command:    1e4*Time[in sec]:
A/=0.5      3.37723994255
A/=0.51     3.42854404449
A*=2.0      1.1587908268
A*2.0       1.19793796539
A+=A        1.11329007149

Теперь нет разницы между делением на .5 и на .51!

Я попробовал его для разных версий numpy и разных машин. На некоторых машинах (например, Intel Xeon E5-2620) этот эффект можно увидеть, но не на некоторых других машинах - и это не зависит от версии numpy.

С помощью script @Ralph Versteegen (см. его отличный ответ!) я получаю следующие результаты:

  • тайминг с i5-2620 (Haswell, 2x6 ядер, но очень старая версия с числовым значением, которая не использует SIMD):

введите описание изображения здесь

  • тайминги с i7-5500U (Broadwell, 2 ядра, numpy 1.11.2):

i7-5500u

Вопрос: В чем причина более высокой стоимости деления на 0.51 по сравнению с делением на 0.5 для некоторых процессоров, если размеры массива велики ( > 10 ^ 6).

В ответе @nneonneo говорится, что для некоторых процессоров Intel существует оптимизация при делении на две степени, но это не объясняет, почему мы можем видеть преимущество этого только для больших массивов.


Первоначальный вопрос: "Как можно объяснить эти разные поведения (деление на 0.5 и деление на 0.51)?"

Здесь также мое первоначальное тестирование script, которое создало тайминги:

import numpy as np
import timeit

def timeit_command( command, rep):
    print "\t"+command+"\t\t", min(timeit.repeat("for i in xrange(%d):"
        %rep+command, "from __main__ import A", number=7))    

sizes=[1e8,  1e4]
reps=[1,  1e4]
commands=["A/=0.5", "A/=0.51", "A*=2.2", "A*=2.0", "A*2.2", "A*2.0",
          "A+=A", "A+A"]

for size, rep in zip(sizes, reps):
    A=np.random.rand(size)
    print "Size:",size
    for command in commands:
        timeit_command(command, rep)

4b9b3361

Ответ 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)

Ответ 2

Процессоры Intel имеют специальную оптимизацию при делении по степеням двух. См., Например, http://www.agner.org/optimize/instruction_tables.pdf, где указано

Задержка FDIV зависит от точности, указанной в контрольном слове: точность 64 бит дает задержку 38, точность 53 бит дает латентность 32, точность 24 бит дает латентность 18. Разделение на мощность 2 занимает 9 часов.

Хотя это применимо к FDIV, а не DIVPD (как примечания ответа @RalphVersteegen), было бы удивительно, если бы DIVPD также не реализовал эту оптимизацию.


Отдел обычно является очень медленным делом. Однако деление на две силы - это просто сдвиг экспоненты, и мантисса обычно не нуждается в изменении. Это делает операцию очень быстрой. Кроме того, легко обнаружить силу двух в представлении с плавающей точкой, так как мантисса будет всеми нулями (с неявным лидирующим 1), поэтому эту оптимизацию легко проверить и дешево реализовать.