为什么在numpy中，log2和log1p比log和log10快得多？

Question

在尝试这个问题时，我注意到了有关np.log2、np.log和np.log10相对性能的差异令我无法解释：

In [1]: %%timeit x = np.random.rand(100000)
   ....: np.log2(x)
   ....: 
1000 loops, best of 3: 1.31 ms per loop

In [2]: %%timeit x = np.random.rand(100000)
np.log(x)
   ....: 
100 loops, best of 3: 3.64 ms per loop

In [3]: %%timeit x = np.random.rand(100000)
np.log10(x)
   ....: 
100 loops, best of 3: 3.93 ms per loop

np.log2比np.log和np.log10快大约3倍。更让人意想不到的是，计算ln(x + 1)的np.log1p(x)与np.log2相当：

In [4]: %%timeit x = np.random.rand(100000)
np.log1p(x)
   ....: 
1000 loops, best of 3: 1.46 ms per loop

我在numpy v1.10.1和v1.8.2中获得了几乎相同的计时结果。

这些运行时性能的差异是否有直观的解释？

- ali_m

---

这个数学SE的答案似乎在说某些方法可以通过“log2”来计算任何对数。这可能意味着np的对数函数的实现在某种程度上依赖于log2和/或ln(x+1)。我认为这与它们的泰勒级数有关。 - R Nar

这是一个非常有趣的观察。我绝不是高效计算例程低级实现方面的专家。直觉上，我猜想这与所有对数在概念上相关有关。如果你知道其中一个，通过简单的转换，你基本上就知道了所有的对数。因此，在某个时候，你必须决定哪一个可以在处理器上高效地计算。通过转换计算其他对数显然需要更多的时间。但我很想看到专家在这里的回答。 - cel

也许由于二进制数据是基于2的，所以有一些log2优化技巧可用。 - wim

这可能与对数函数log(x+1)的Taylor级数相对简单有关。 - R Nar

@FermionPortal，您有兴趣将您的评论写成答案吗？我可以自己试一下，但让赏金白白浪费似乎有点可惜;-) - ali_m

Answer 1

我发现您的问题非常有趣，所以花了几个小时进行了一些研究; 我认为我已经找到了一个关于性能差异的解释适用于整数（感谢Matteo Italia的提醒）- 尽管如此，不清楚这种推理如何扩展到浮点数：

计算机使用基数2-根据参考链接中的文章，计算log2是一个需要4个处理器周期的过程-计算log10需要将log2（val）乘以1 / log2（10），这又需要5个周期。

找到log2是查找值的最低有效位的索引的问题。 （视频大约在第23分钟开始）。

位运算技巧：找到整数的以10为底的对数

位操作技巧：在O(lg(N))时间内找到N位整数的对数基数

整数以10为底的对数是通过先使用上述方法之一来找到以2为底的对数。由于关系log10(v) = log2(v) / log2(10)，我们需要将其乘以1/log2(10)，约为1233/4096，或者是1233后面跟着右移12位。加1是必要的，因为IntegerLogBase2向下取整。最后，由于值t只是可能偏差1的近似值，所以准确值是通过减去v < PowersOf10[t]的结果获得的。

这种方法比IntegerLogBase2多6个操作。它可以通过修改上面的以2为底的查找表格方法，使条目保持t的计算结果（即预添加、-乘法和-移位）来加速（在具有快速内存访问的机器上）。这样做只需要总共9个操作来找到以10为底的对数，假设使用了4个表格（每个字节一个表格）。

注意：在MIT视频：软件系统性能工程，2010年秋季|Lec 2 | MIT 6.172 Performance Engineering of Software Systems, Fall 2010（从第36分钟开始）中使用了DeBruijn序列查找和位移技术来计算log2。

请注意，这篇StackOverflow帖子演示了一种使高效的C++ log2计算真正跨平台的方法。

警告：我没有检查numpy源代码以验证它是否实现了类似的技术，但如果它没有实现，那将令人惊讶。事实上，从OP帖子下面的评论中，Fermion Portal已经进行了检查：

实际上，numpy使用glibc中的math.h，如果您使用math.h/cmath.h在C/C++中，您将看到相同的差异。您可以找到两个函数的注释良好的源代码，例如ic.unicamp.br/~islene/2s2008-mo806/libc/sysdeps/ieee754/dbl-64/…和ic.unicamp.br/~islene/2s2008-mo806/libc/sysdeps/ieee754/dbl-64/… - Fermion Portal[9]

- Reblochon Masque

---

注意，这里我们谈论的是浮点数上的对数运算，上面的“位操作技巧”不适用。 - Matteo Italia

哦糟糕！啊！就在我认为我的贡献可能有用的时候！O_O...我会添加一条注释，说明它适用于整数；也许对某些人仍然有用？谢谢你指出来，马特奥 - 我还是学到了东西！ :) - Reblochon Masque

Answer 2

这只是一条注释，但比评论更长。显然与您的特定安装有关：

import numpy as np
import numexpr as ne
x = np.random.rand(100000)

我使用conda中的numpy 1.10和使用icc编译的版本得到了相同的时间：

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

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

我认为这可能与获取MKL VML软件包有关，但看起来并不是:

%timeit ne.evaluate('log(x)')
1000 loops, best of 3: 218 µs per loop

看起来你的numpy安装从两个不同的地方获取其log / log2实现，这很奇怪。

- Daniel

---

很有趣-我仍然在使用numpy 1.10.1或conda中的1.9.3看到np.log和np.log2非常不同的时间，尽管这两个似乎都是使用gcc 4.4.7编译而成，而不是icc。我无法访问使用icc编译的版本以进一步测试。 - ali_m

我获得了 ICC 16.0.1 的副本，并按照此处的说明从源代码构建了 numpy 1.10.1。现在我发现整体性能稍微变差了，使用 np.log 和 np.log10 仍然比使用 np.log2 和 np.log1p 慢大约 2 倍。 - ali_m

@ali_m 更加好奇。你是否碰巧在使用 AMD CPU？ - Daniel

没有 - 到目前为止，我已经在两台英特尔机器上尝试过了。你的设置是什么？ - ali_m

@ali_m 也是英特尔机器。你试过使用 conda 吗？ - Daniel

是的 - 正如我上面提到的，我也尝试过通过 conda 安装的 numpy 1.10.1 和 1.9.3。我还尝试了它们与可选的 MKL 优化一起使用和不使用。在每种情况下，我都看到 np.log2 和 np.log 之间至少有两倍的差异。 - ali_m

Answer 3

免责声明：我既不是可靠的，也不是官方的来源。

我几乎可以确定，任何以自然对数为底数的对数函数的实现都可以与log2函数一样快，因为要将一个转换为另一个，只需要除以一个常数。当然，这是在假设单个除法操作只占其他计算的一小部分的情况下成立的；在对数的准确实现中，这是正确的。

在大多数情况下，numpy使用glibc的math.h，如果你使用math.h/cmath.h，你会在C/C++中看到相同的差异。在评论中，一些人观察到np.log和np.log2具有相同的速度；我怀疑这可能来自不同的构建/平台。

你可以在this GitHub repo的dbl-64/和flt-32/子目录中的e_log.c、e_log2.c、e_logf.c和e_log2f.c文件中找到这两个函数的注释良好的源代码。

对于双精度，在glibc中，log函数实现了一个与IBM从2001年左右的libultim库中不同的算法（与log2相比），而log2来自Sun Microsystems约1993年。仅仅看一下代码，就可以看出正在实现两种不同的近似值。相比之下，对于单精度，log和log2这两个函数除了在log2情况下除以ln2之外完全相同，因此速度相同。

关于底层算法、替代方案以及讨论哪些算法应该在未来的glibc中包含，请参见此处。

- Fermion Portal

---

感谢您深入研究此问题。我将授予您赏金，因为我认为这些链接迄今为止提供了最有用的洞察力。然而，我仍然犹豫是否将此问题标记为关闭，以防其他人想要提供更全面的答案。 - ali_m

Answer 4

这可能应该是一个注释，但会太长...

为了使内容更有趣，在2018年的Windows 10 64位机器上，结果被颠倒了。

默认Anaconda

Python 3.6.3 |Anaconda, Inc.| (default, Oct 15 2017, 03:27:45) [MSC v.1900 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 6.1.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np; np.random.seed(0); x = np.random.rand(100000)
   ...: %timeit np.log2(x)
   ...: %timeit np.log1p(x)
   ...: %timeit np.log(x)
   ...: %timeit np.log10(x)
   ...:
1.48 ms ± 18 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
1.33 ms ± 36.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
840 µs ± 7.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
894 µs ± 2.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Intel Python

Python 3.6.3 |Intel Corporation| (default, Oct 17 2017, 23:26:12) [MSC v.1900 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 6.1.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np; np.random.seed(0); x = np.random.rand(100000)
   ...: %timeit np.log2(x)
   ...: %timeit np.log1p(x)
   ...: %timeit np.log(x)
   ...: %timeit np.log10(x)
   ...:
1.01 ms ± 2.57 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
236 µs ± 6.08 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
161 µs ± 1.77 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
171 µs ± 1.34 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

- Constructor