与 astype(int) 相比，numpy around/rint 慢答案

【问题标题】：numpy around/rint slow compared to astype(int)与 astype(int) 相比，numpy around/rint 慢
【发布时间】：2015-01-30 20:44:34
【问题描述】：

所以如果我有类似x=np.random.rand(60000)*400-200 的东西。 iPython 的 %timeit 说：

x.astype(int) 耗时 0.14 毫秒
np.rint(x) 和 np.around(x) 需要 1.01 毫秒

请注意，在 rint 和 around 情况下，您仍然需要花费额外的 0.14 毫秒来完成最终的 astype(int)（假设这是您最终想要的）。

问题：我认为大多数现代硬件能够同时执行这两项操作是否正确。如果是这样，为什么 numpy 的舍入时间要长 8 倍？

碰巧我对算术的准确性并不是特别挑剔，但我不知道如何利用 numpy 来利用这一点（我正在做的是混乱的生物学而不是粒子物理学）。

【问题讨论】：

承认我对 CPU 的算术能力不是很熟悉：为什么他们能够在相同的时间内完成它？ astype 只是切除了一些位，舍入操作需要检查你切除了多少（以确定你是否舍入到较低或较高的 int）。
嗯，有一件事是从浮点数转换为整数类型只需丢弃小数部分，这相当于向零舍入，而np.rint 则舍入到最近的 整数（这是额外的工作）。因此，np.trunc 与np.astype(int) 更具可比性。在我的速度测试中，np.trunc 仍然较慢，但是查看source，这可能是因为它是根据 ceil 和 floor 实现的，而不是简单的演员表。
似乎有低级标志来控制舍入模式，例如：gcc.gnu.org/wiki/FloatingPointMath。

标签： python c assembly numpy sse

【解决方案1】：

np.around(x).astype(int) 和 x.astype(int) 不会产生相同的值。前轮偶数（与((x*x>=0+0.5) + (x*x<0-0.5)).astype(int) 相同），而后轮趋向于零。然而，

y = np.trunc(x).astype(int)
z = x.astype(int)

显示y==z，但计算y 要慢得多。所以是 np.trunc 和 np.around 函数很慢。

In [165]: x.dtype
Out[165]: dtype('float64')
In [168]: y.dtype
Out[168]: dtype('int64')

所以np.trunc(x) 从两倍到两倍向零舍入。然后astype(int) 必须将 double 转换为 int64。

在内部，我不知道 python 或 numpy 在做什么，但我知道如何在 C 中做到这一点。让我们讨论一些硬件。使用 SSE4.1，可以使用以下方法从 double 到 double 进行 round、floor、ceil 和 trunc：

_mm_round_pd(a, 0); //round: round even
_mm_round_pd(a, 1); //floor: round towards minus infinity
_mm_round_pd(a, 2); //ceil:  round towards positive infinity
_mm_round_pd(a, 3); //trunc: round towards zero

但 numpy 也需要支持没有 SSE4.1 的系统，因此它必须在没有 SSE4.1 和 SSE4.1 的情况下构建，然后使用调度程序。

但是使用 SSE/AVX 从 double 直接到 int64 执行此操作直到 AVX512 才有效。但是，可以仅使用 SSE2 有效地将 double 舍入为 int32：

_mm_cvtpd_epi32(a);  //round double to int32 then expand to int64
_mm_cvttpd_epi32(a); //trunc double to int32 then expand to int64

这些将两个双精度数转换为两个 int64。

在您的情况下，这可以正常工作，因为范围肯定在 int32 内。但是除非 python 知道范围适合 int32，否则它不能假设这一点，所以它必须舍入或截断到 int64，这很慢。此外，无论如何，必须再次构建 numpy 以支持 SSE2。

但也许您可以使用单个浮点数组开始。在这种情况下，您可以这样做：

_mm_cvtps_epi32(a); //round single to int32
_mm_cvttps_epi32(a) //trunc single to int32

这些将四个单曲转换为四个 int32。

因此，要回答您的问题，SSE2 可以有效地从 double 舍入或截断为 int32。 AVX512 将能够使用_mm512_cvtpd_epi64(a) 或_mm512_cvttpd_epi64(a) 有效地将双精度整数舍入或截断为int64。 SSE4.1可以高效地从float到float或double到double的round/trunc/floor/ceil。

【讨论】：

感谢您提供详细信息。我刚刚尝试将x 设为单一并使用int32 而不是int...但没有运气。我还注意到rint 有一个out arg，您可以将其设置为一个空的int32 数组，但这对速度也没有帮助。（我无法确定我的 numpy 是用哪个 SSE 版本编译的，但它可能是 SSE2 或 SSE3，因为它是一台足够好的机器，而且 Windows 安装程序应该没有 SSE、SSE2 和 SSE3 二进制文件。）
@dan-man，是的，我尝试了 np.float32 和 np.int32 以及其他变体，但没有任何改进。但这对你来说真的是个问题吗？ “过早的优化是万恶之源”。这是 numpy 开发人员应该担心的事情。
我希望一个 numpy 开发人员会出现并告诉我一个快速的 hack 或指出一个已知的错误......如果这样的话那将是值得的，因为我有一个花费 1 秒的函数（ >50% 总时间）在rint 上，我需要调用该函数数百次。但是，是的，我可能应该承认某种程度的过早优化。（虽然，当我发现这种差异时，我想到可能有一整类这样的事情，其中 numpy 比它应该慢一个数量级......如果是这样，那么值得了解它们。）
@dan-man，在这种情况下，您可能希望在 SO 上发布您的功能并查看您得到的答案。我通常在 C 中进行优化。我们中的一些人在工作中也使用 pyopencl。您可能要考虑这一点。我过去使用过它，但现在使用 C 语言在 CPU 上使用 OpenMP 和 SIMD 可以获得足够好的结果。

【解决方案2】：

正如@jme 在 cmets 中指出的那样，rint 和 around 函数必须确定是向上还是向下舍入到最接近的整数。相反，astype 函数将始终向下舍入，因此它可以立即丢弃十进制信息。还有许多其他功能可以做同样的事情。此外，您可以通过使用较少的整数位数来提高速度。但是，您必须小心，您可以容纳输入数据的全部范围。

%%timeit
np.int8(x)
10000 loops, best of 3: 165 µs per loop

注意，这不会存储 -128 到 127 范围之外的值，因为它是 8 位的。您的示例中的某些值超出了此范围。

在我尝试过的所有其他方法中，np.intc 似乎是最快的：

%%timeit
np.int16(x)
10000 loops, best of 3: 186 µs per loop

%%timeit
np.intc(x)
10000 loops, best of 3: 169 µs per loop

%%timeit
np.int0(x)
10000 loops, best of 3: 170 µs per loop

%%timeit
np.int_(x)
10000 loops, best of 3: 188 µs per loop

%%timeit
np.int32(x)
10000 loops, best of 3: 187 µs per loop

%%timeit
    np.trunc(x)
1000 loops, best of 3: 940 µs per loop

你的例子，在我的机器上：

%%timeit
np.around(x)
1000 loops, best of 3: 1.48 ms per loop

%%timeit
np.rint(x)
1000 loops, best of 3: 1.49 ms per loop

%%timeit
x.astype(int)
10000 loops, best of 3: 188 µs per loop

【讨论】：