正弦计算比余弦慢几个数量级

Posted 2023-02-14

【中文标题】正弦计算比余弦慢几个数量级

【英文标题】：sine calculation orders of magnitude slower than cosine

【发布时间】：2016-06-19 07:29:42

【问题描述】：

tl;博士

对于相同的numpy 数组，计算np.cos 需要3.2 秒，而np.sin 在Linux Mint 上运行548 秒（九分钟）。

完整代码见this repo。

---

我有一个脉冲信号（见下图），我需要将其调制到 HF 载波上，模拟Laser Doppler Vibrometer。因此需要对信号及其时基进行重采样以匹配载波更高的采样率。

在接下来的解调过程中，同相载波cos(omega * t)和相移载波sin(omega * t)都是需要的。 奇怪的是，评估这些函数的时间很大程度上取决于计算时间向量的方式。

时间向量t1直接使用np.linspace计算，t2使用method implemented in 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]

如下图所示，时间向量并不相同。在 8000 万个样本中，t1 - t2 的差异达到了 1e-8。

在我的机器上计算t1 的同相和偏移载波分别需要 3.2 秒。但是，对于 t2，计算偏移载波需要 >540 秒。九分钟。对于几乎相同的 8000 万个值。

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

我可以在运行 Linux Mint 17 的 32 位笔记本电脑和 64 位塔式电脑上重现此错误。然而，在我室友的 MacBook 上，“慢正弦”所需的时间与其他三个计算一样少。

---

我在 64 位 AMD 处理器上运行 Linux Mint 17.03，在 32 位 Intel 处理器上运行 Linux Mint 17.2。

【问题讨论】：

如果切换调用这些的顺序，结果是否一致？ （只是假设这可能是由于某种内存/缓存问题 - 每一个都会产生一个 640MB 的向量。）

@OliverCharlesworth 是的，结果是一样的

numpy 链接的库是否在计算机之间有所不同？

@MSeifert 我该如何检查？

numpy.__config__.show()

【参考方案1】：

造成这些巨大性能差异的一个可能原因可能是数学库如何创建或处理 IEEE 浮点下溢（或去范数），这可能是由超越函数逼近过程中一些较小尾数位的差异产生的。并且您的 t1 和 t2 向量可能会因这些较小的尾数位以及您链接的任何库中用于计算超越函数的算法以及每个特定操作系统上的 IEEE 算术 denorms 或下溢处理程序而有所不同。

【参考方案2】：

我认为 numpy 与此无关：我认为您在系统上的 C 数学库中遇到了性能错误，该错误会影响 pi 的大倍数附近的 sin。 （我在这里使用了广义上的“bug”——据我所知，由于大浮点数的正弦定义不明确，“bug”实际上是库在处理极端情况时表现正确！）

在 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

以及Python chatroom 报告中的其他使用 Linux 的类型

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

【讨论】：

这可能与denormal numbers 有关吗？我记得写过一些浮点代码，当它非常小但涉及非零数字时会变得非常慢。

但是这个“错误”对 大 数字生效，而不是接近于零的数字...

@Finwood 不是解释，但如果问题是大数字，你能把它当作 mod 2pi 吗？

@Paul 是的，这就是我正在做的规避问题