27

在 NumPy 中，x*x*x 比 x**3 甚至 np.power(x, 3) 快一个数量级。

x = np.random.rand(1e6)
%timeit x**3
100 loops, best of 3: 7.07 ms per loop

%timeit x*x*x
10000 loops, best of 3: 163 µs per loop

%timeit np.power(x, 3)
100 loops, best of 3: 7.15 ms per loop

关于为什么会发生这种行为的任何想法？据我所知，这三个产生相同的输出（用 np.allclose 检查）。

4

6 回答 6

36

根据这个答案，这是因为求幂的实现有一些乘法没有的开销。然而，随着指数的增加，朴素的乘法会变得越来越慢。经验证明：

 In [3]: x = np.random.rand(1e6)

 In [15]: %timeit x**2
 100 loops, best of 3: 11.9 ms per loop

 In [16]: %timeit x*x
 100 loops, best of 3: 12.7 ms per loop

 In [17]: %timeit x**3
 10 loops, best of 3: 132 ms per loop

 In [18]: %timeit x*x*x
 10 loops, best of 3: 27.2 ms per loop

 In [19]: %timeit x**4
 10 loops, best of 3: 132 ms per loop

 In [20]: %timeit x*x*x*x
 10 loops, best of 3: 42.4 ms per loop

 In [21]: %timeit x**10
 10 loops, best of 3: 132 ms per loop

 In [22]: %timeit x*x*x*x*x*x*x*x*x*x
 10 loops, best of 3: 137 ms per loop

 In [24]: %timeit x**15
 10 loops, best of 3: 132 ms per loop

 In [25]: %timeit x*x*x*x*x*x*x*x*x*x*x*x*x*x*x
 1 loops, best of 3: 212 ms per loop

请注意，求幂时间或多或少保持不变，除了x**2我怀疑是特殊情况的情况，而乘法变得越来越慢。看来您可以利用它来获得更快的整数幂运算......例如：

In [26]: %timeit x**16
10 loops, best of 3: 132 ms per loop

In [27]: %timeit x*x*x*x*x*x*x*x*x*x*x*x*x*x*x*x
1 loops, best of 3: 225 ms per loop

In [28]: def tosixteenth(x):
   ....:     x2 = x*x
   ....:     x4 = x2*x2
   ....:     x8 = x4*x4
   ....:     x16 = x8*x8
   ....:     return x16
   ....:

In [29]: %timeit tosixteenth(x)
10 loops, best of 3: 49.5 ms per loop

看来您可以通过将任何整数拆分为 2 的幂和，计算上述 2 的每个幂，然后求和来一般地应用此技术：

In [93]: %paste
def smartintexp(x, exp):
    result = np.ones(len(x))
    curexp = np.array(x)
    while True:
        if exp%2 == 1:
            result *= curexp
        exp >>= 1
        if not exp: break
        curexp *= curexp
    return result
## -- End pasted text --

In [94]: x
Out[94]:
array([ 0.0163407 ,  0.57694587,  0.47336487, ...,  0.70255032,
        0.62043303,  0.0796748 ])

In [99]: x**21
Out[99]:
array([  3.01080670e-38,   9.63466181e-06,   1.51048544e-07, ...,
         6.02873388e-04,   4.43193256e-05,   8.46721060e-24])

In [100]: smartintexp(x, 21)
Out[100]:
array([  3.01080670e-38,   9.63466181e-06,   1.51048544e-07, ...,
         6.02873388e-04,   4.43193256e-05,   8.46721060e-24])

In [101]: %timeit x**21
10 loops, best of 3: 132 ms per loop

In [102]: %timeit smartintexp(x, 21)
10 loops, best of 3: 70.7 ms per loop

对于两个的小偶数幂，它很快：

In [106]: %timeit x**32
10 loops, best of 3: 131 ms per loop

In [107]: %timeit smartintexp(x, 32)
10 loops, best of 3: 57.4 ms per loop

但随着指数变大变慢：

In [97]: %timeit x**63
10 loops, best of 3: 133 ms per loop

In [98]: %timeit smartintexp(x, 63)
10 loops, best of 3: 110 ms per loop

对于大的最坏情况并没有更快：

In [115]: %timeit x**511
10 loops, best of 3: 135 ms per loop

In [114]: %timeit smartintexp(x, 511)
10 loops, best of 3: 192 ms per loop
于 2013-08-26T22:21:46.277 回答
7

如果您正在计算功率并且担心速度，请注意：

x = np.random.rand(5e7)

%timeit x*x*x
1 loops, best of 3: 522 ms per loop

%timeit np.einsum('i,i,i->i',x,x,x)
1 loops, best of 3: 288 ms per loop

为什么 einsum 更快仍然是我的一个悬而未决的问题。虽然它喜欢由于einsum能够使用 SSE2 而 numpy 的 ufuncs 直到 1.8 才会使用。

就地更快：

def calc_power(arr):
    for x in xrange(arr.shape[0]):
        arr[x]=arr[x]*arr[x]*arr[x]
numba_power = autojit(calc_power)

%timeit numba_power(x)
10 loops, best of 3: 51.5 ms per loop

%timeit np.einsum('i,i,i->i',x,x,x,out=x)
10 loops, best of 3: 111 ms per loop

%timeit np.power(x,3,out=x)
1 loops, best of 3: 609 ms per loop
于 2013-08-26T22:31:29.430 回答
3

我希望这是因为必须处理两者都是浮动x**y的通用情况。数学上我们可以写. 按照你的例子，我发现xyx**y = exp(y*log(x))

x = np.random.rand(1e6)
%timeit x**3
10 loops, best of 3: 178 ms per loop

%timeit np.exp(3*np.log(x))
10 loops, best of 3: 176 ms per loop

我没有检查实际的 numpy 代码，但它必须在内部执行类似的操作。

于 2013-08-26T22:19:46.387 回答
1

这是因为 python 中的幂是作为浮点运算执行的（对于 numpy 也是如此，因为它使用 C）。

在 C 中，pow 函数提供了 3 种方法：

双 pow（双 x，双 y）

长碗（长双 x，长双 y）

浮动 powf（浮动 x，浮动 y）

其中每一个都是浮点运算。

于 2013-08-26T22:12:10.447 回答
0

根据规范：

两个参数的形式 pow(x, y) 等价于使用幂运算符：x**y。

参数必须具有数字类型。对于混合操作数类型，适用二元算术运算符的强制规则。

换句话说：因为x是浮点数，所以指数从整数转换为浮点数，并执行通用浮点幂运算。在内部，这通常被重写为：

x**y = 2**(y*lg(x))

2**a和lg a（以 2 为底的对数a）是现代处理器上的单条指令，但它仍然需要比几次乘法更长的时间。

于 2013-08-26T22:36:04.677 回答
-1 
timeit np.multiply(np.multiply(x,x),x)

次相同x*x*x。我的猜测是np.multiply使用像 BLAS 这样的快速 Fortran 线性代数包。我从另一个numpy.dot在某些情况下使用 BLAS 的问题中知道。

我必须把它拿回来。 np.dot(x,x)比 快 3 倍np.sum(x*x)。所以速度优势np.multiply与使用 BLAS 不一致。

使用我的 numpy（时间会因机器和可用库而异）

np.power(x,3.1)
np.exp(3.1*np.log(x))

大约需要相同的时间，但

np.power(x,3)

是 2 倍的速度。速度不如x*x*x，但还是比一般的功率快。所以它利用了整数幂的一些优势。

于 2013-08-26T22:29:46.540 回答