numpy.cos работает значительно дольше на определенных номерах

TL;DR:

numpy.cos() работает на 30% дольше для определенных чисел (например, ровно 24000.0). Добавление небольшой дельты (+0.01) приводит к numpy.cos() работе numpy.cos().

Понятия не имею почему.

Я наткнулся на странную проблему во время моей работы с numpy. Я проверял работу кеша и случайно сделал неправильный график - как время numpy.cos(X) зависит от X Вот мой модифицированный код (скопированный из моего блокнота Jupyter):

import numpy as np
import timeit
st = 'import numpy as np'
cmp = []
cmp_list = []
left = 0
right = 50000
step = 1000
# Loop for additional average smoothing
for _ in range(10):
    cmp_list = []
    # Calculate np.cos depending on its argument
    for i in range(left, right, step):
        s=(timeit.timeit('np.cos({})'.format(i), number=15000, setup=st))
        cmp_list.append(int(s*1000)/1000)
    cmp.append(cmp_list)

# Calculate average times
av=[np.average([cmp[i][j] for i in range(len(cmp))]) for j in range(len(cmp[0]))]

# Draw the graph
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.plot(range(left, right, step), av, marker='.')
plt.show()

График выглядел так:

enter image description here

Сначала я подумал, что это просто случайный сбой. Я пересчитал свои камеры, но результат был почти таким же. Поэтому я начал играть с параметрами step, с количеством вычислений и средней длиной списка. Но все это никак не повлияло на это число:

enter image description here

И еще ближе:

enter image description here

После этого range был бесполезен (он не может np.cos с поплавками), поэтому я вычислил np.cos вручную:

print(timeit.timeit('np.cos({})'.format(24000.01),number=5000000,setup=st))
print(timeit.timeit('np.cos({})'.format(24000.00),number=5000000,setup=st))
print(timeit.timeit('np.cos({})'.format(23999.99),number=5000000,setup=st))

И результат был:

3.4297256958670914
4.337243619374931
3.4064380447380245

np.cos() рассчитывает ровно 24000,00 на 30% дольше, чем 24000.01!

Были и другие странные цифры вроде этого (где-то около 500000, я точно не помню).

Я просмотрел numpy документацию, ее исходный код, и он ничего не имел об этом эффекте. Я знаю, что тригонометрические функции используют несколько алгоритмов, в зависимости от размера значения, точности и т.д., Но меня смущает, что точные числа можно вычислить намного дольше.

Почему np.cos() имеет такой странный эффект? Это какой-то побочный эффект процессора (потому что numpy.cos использует C-функции, которые зависят от процессоров)? У меня установлены Intel Core i5 и Ubuntu, если это кому-то поможет.

Редактировать 1: я пытался воспроизвести его на другом компьютере с AMD Ryzen 5. Результаты просто нестабильны. Вот графики для трех последовательных прогонов одного и того же кода:

import numpy as np
import timeit

s = 'import numpy as np'
times = []
x_ranges = np.arange(23999, 24001, 0.01)
for x in x_ranges:
    times.append(timeit.timeit('np.cos({})'.format(x), number=100000, setup=s))

# ---------------

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(x_ranges, times)
plt.show()

enter image description here

enter image description here

enter image description here

Ну, есть некоторые паттерны (например, в основном непротиворечивая левая часть и непоследовательная правая часть), но она сильно отличается от процессоров Intel. Похоже, что это действительно особые аспекты процессоров, а поведение AMD гораздо более предсказуемо в своей неопределенности :)

PS @WarrenWeckesser спасибо за функцию '' np.arange '' '. Это действительно полезно, но это ничего не меняет в результатах, как и ожидалось.

Ответы

Ответ 1

Медлительность вычисления результата для этих специальных чисел может быть связана с точной дилеммой округления и создания таблицы.

Для иллюстрации предположим, что вы составляете таблицу экспоненциальной функции на 4 места. Тогда exp (1,626) = 5,0835. Должно ли это быть округлено до 5,083 или 5,084? Если exp (1.626) вычисляется более тщательно, он становится 5.08350. А потом 5.083500. А потом 5.0835000. Поскольку exp трансцендентен, это может продолжаться произвольно задолго до того, как будет определено, является ли exp (1.626) 5.083500... 0ddd или 5.0834999... 9ddd.

Хотя по этой причине стандарт IEEE НЕ требует, чтобы трансцендентные функции были точно округлены, возможно, что реализация функции math.cos этой проблемой, math.cos все усилия для вычисления наиболее точного результата, а затем выяснения что эффект не стоит усилий.

Чтобы продемонстрировать, что это так для некоторого числа X, нужно будет с высокой точностью вычислить значение math.cos(X) и проверить его двоичное представление - за представительной частью мантиссы должен следовать один из следующих шаблонов:

  • 1 и длинный пробег 0
  • 0 и длительный прогон 1 (когда значение вычисляется с точностью ниже, чем требуется для размещения всех 1 в прогоне, этот случай отображается как первый)

Следовательно, вероятность того, что число будет медленным аргументом для трансцендентной функции, равна 1/2 n, где n - максимальная длина вышеуказанного паттерна, видимого алгоритмом, после которого он прекращает попытки получить точно округленный результат.

Демонстрация с выделением представительной части мантиссы для случая двойной точности IEEE 754 (где у мантиссы 53 бита):

In [1]: from mpmath import mp

In [2]: import math

In [3]: def show_mantissa_bits(x, n, k):
   ...:     print(bin(int(mp.floor(abs(x) * 2**n)))[2:])
   ...:     print('^'*k)
   ...:     

In [4]: mp.prec = 100

In [5]: show_mantissa_bits(mp.cos(108), 64, 53)
110000000100001011001011010000110111110010100011000011000000000
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [6]: show_mantissa_bits(mp.cos(108.01), 64, 53)
101110111000000110001101110001000010100111000010101100000100110
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [7]: show_mantissa_bits(mp.cos(448), 64, 53)
101000101000100111000010111100001011111000001111110001000000000
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [8]: show_mantissa_bits(mp.cos(448.01), 64, 53)
101001110110001010010100100000110001111100000001101110111010111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [9]: show_mantissa_bits(mp.cos(495), 64, 53)
11001010100101110110001100110101010011110010000000000011111111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [10]: show_mantissa_bits(mp.cos(495.01), 64, 53)
11010100100111100110000000011000110000001001101100010000001010
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [11]: show_mantissa_bits(mp.cos(24000), 64, 53)
11001000100000001100110111011101001101101101000000110011111111
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In [12]: show_mantissa_bits(mp.cos(24000.01), 64, 53)
10111110011100111001010101100101110001011010101011001010110011
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^