Вычисление порядка синуса медленнее, чем косинус

TL;DR

Из того же массива numpy, вычисление np.cos занимает 3,2 секунды, если np.sin работает 548 секунд (девять минут) в Linux Mint.

См. это репо для полного кода.

У меня есть импульсный сигнал (см. изображение ниже), который мне нужно модулировать на HF-носитель, имитируя Лазерный допплеровский виброметр. Поэтому сигнал и его временной интервал должны быть повторно сэмплированы в соответствии с более высокой частотой дискретизации несущей.

В следующем процессе демодуляции необходимы как фазовая несущая cos(omega * t), так и несущая фаза sin(omega * t). Как ни странно, время для оценки этих функций сильно зависит от того, как был рассчитан вектор времени.

. Временной вектор t1 вычисляется с использованием np.linspace напрямую, t2 использует метод реализованный в scipy.signal.resample.

pulse = np.load('data/pulse.npy')  # 768 samples

pulse_samples = len(pulse)
pulse_samplerate = 960  # 960 Hz
pulse_duration = pulse_samples / pulse_samplerate  # here: 0.8 s
pulse_time = np.linspace(0, pulse_duration, pulse_samples,
                         endpoint=False)

carrier_freq = 40e6  # 40 MHz
carrier_samplerate = 100e6  # 100 MHz
carrier_samples = pulse_duration * carrier_samplerate  # 80 million

t1 = np.linspace(0, pulse_duration, carrier_samples)

# method used in scipy.signal.resample
# https://github.com/scipy/scipy/blob/v0.17.0/scipy/signal/signaltools.py#L1754
t2 = np.arange(0, carrier_samples) * (pulse_time[1] - pulse_time[0]) \
        * pulse_samples / float(carrier_samples) + pulse_time[0]

Как видно на рисунке ниже, векторы времени не идентичны. При 80 миллионах выборок разность t1 - t2 достигает 1e-8.

Вычисление синфазной и сдвинутой несущей t1 занимает на каждой машине 3,2 секунды.
С t2, однако, вычисление сдвинутой несущей занимает 540 секунд. Девять минут. Почти для тех же 80 миллионов значений.

omega_t1 = 2 * np.pi * carrier_frequency * t1
np.cos(omega_t1)  # 3.2 seconds
np.sin(omega_t1)  # 3.3 seconds

omega_t2 = 2 * np.pi * carrier_frequency * t2
np.cos(omega_t2)  # 3.2 seconds
np.sin(omega_t2)  # 9 minutes

Я могу воспроизвести эту ошибку как на моем 32-битном ноутбуке, так и на моей 64-битной башне, работающей под управлением Linux Mint 17. Однако на моем плоском матке MacBook "медленный синус" занимает столько же времени, сколько и остальные три вычисления.

Я запускаю Linux Mint 17.03 на 64-разрядном процессоре AMD и Linux Mint 17.2 на 32-разрядном процессоре Intel.

Ответ 1

Я не думаю, что numpy имеет какое-то отношение к этому: я думаю, что вы отключите ошибку производительности в C-математической библиотеке вашей системы, которая влияет на грех вблизи больших кратных значений pi. (Я использую "ошибку" в довольно широком смысле здесь - насколько я знаю, поскольку синус больших поплавков плохо определен, "ошибка" на самом деле является библиотекой, которая ведет себя правильно, чтобы обрабатывать угловые случаи!)

В linux я получаю:

>>> %timeit -n 10000 math.sin(6e7*math.pi)
10000 loops, best of 3: 191 µs per loop
>>> %timeit -n 10000 math.sin(6e7*math.pi+0.12)
10000 loops, best of 3: 428 ns per loop

и другие типы использования Linux из отчета Python:

10000 loops, best of 3: 49.4 µs per loop 
10000 loops, best of 3: 206 ns per loop

In [3]: %timeit -n 10000 math.sin(6e7*math.pi)
10000 loops, best of 3: 116 µs per loop

In [4]: %timeit -n 10000 math.sin(6e7*math.pi+0.12)
10000 loops, best of 3: 428 ns per loop

но пользователь Mac сообщил

In [3]: timeit -n 10000 math.sin(6e7*math.pi)
10000 loops, best of 3: 300 ns per loop

In [4]: %timeit -n 10000 math.sin(6e7*math.pi+0.12)
10000 loops, best of 3: 361 ns per loop

без разницы по порядку величины. В качестве обходного пути вы можете сначала попробовать вещи mod 2 pi:

>>> new = np.sin(omega_t2[-1000:] % (2*np.pi))
>>> old = np.sin(omega_t2[-1000:])
>>> abs(new - old).max()
7.83773902468434e-09

который имеет лучшую производительность:

>>> %timeit -n 1000 new = np.sin(omega_t2[-1000:] % (2*np.pi))
1000 loops, best of 3: 63.8 µs per loop
>>> %timeit -n 1000 old = np.sin(omega_t2[-1000:])
1000 loops, best of 3: 6.82 ms per loop

Обратите внимание, что, как и ожидалось, аналогичный эффект происходит для cos, только сдвинутый:

>>> %timeit -n 1000 np.cos(6e7*np.pi + np.pi/2)
1000 loops, best of 3: 37.6 µs per loop
>>> %timeit -n 1000 np.cos(6e7*np.pi + np.pi/2 + 0.12)
1000 loops, best of 3: 2.46 µs per loop

Ответ 2

Одна из возможных причин этих огромных различий в производительности может заключаться в том, как математическая библиотека создает или обрабатывает нисходящую точку с плавающей запятой IEEE (или denorms), которая может быть вызвана различием некоторых из более мелких битов мантиссы во время приближения трансцендентальной функции. И ваши векторы t1 и t2 могут отличаться этими меньшими битами мантиссы, а также алгоритм, используемый для вычисления трансцендентальной функции в любых связанных с вами библиотеках, а также арифметические дескрипторы IEEE или обработчик нижнего потока на каждой конкретной ОС.