向量化Python中重复的数学函数

Question

我有一个形式为$f(x)=\sum_{j=0}^N x^j * \sin(j*x)$的数学函数，我想用Python高效地计算它。N是~100阶的。对于一个巨大矩阵的所有条目x，这个函数f都会被计算数千次，因此我希望提高性能(profiler表明f的计算占用了大部分时间)。为了避免在函数f的定义中出现循环，我写道：

def f(x)
    J=np.arange(0,N+1)
    return sum(x**J*np.sin(j*x))

问题是，如果我想为矩阵的所有元素计算这个函数，我需要首先使用numpy.vectorize，但据我所知，这不一定比for循环快。

有没有一种有效的方法来执行这种类型的计算？

Answer 1

欢迎来到Sack Overflow！^^

嗯，计算something ** 100是一件很严肃的事情。但请注意，当您声明数组J时，您是如何强制您的函数独立计算x, x^2, x^3, x^4, ... (依此类推)的。

让我们以这个函数为例(这就是您正在使用的函数)：

def powervector(x, n):
    return x ** np.arange(0, n)

现在还有另一个函数，它甚至不使用NumPy：

def power(x, n):
    result = [1., x]
    aux = x
    for i in range(2, n):
        aux *= x               
        result.append(aux)
    return result

现在，让我们验证它们是否都计算相同的东西：

In []: sum(powervector(1.1, 10))
Out[]: 15.937424601000005

In []: sum(power(1.1, 10))
Out[]: 15.937424601000009

酷，现在让我们比较一下两者的性能(在iPython中)：

In [36]: %timeit sum(powervector(1.1, 10))
The slowest run took 20.42 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 3.52 µs per loop

In [37]: %timeit sum(power(1.1, 10))
The slowest run took 5.28 times longer than the fastest. This could mean that an intermediate result is being cached 
1000000 loops, best of 3: 1.13 µs per loop

它更快，因为你没有计算x的所有能力，因为你知道x ^ N == (x ^ N - 1) * x和你利用它。

您可以使用它来查看您的性能是否有所提高。当然，您可以更改power()以使用NumPy向量作为输出。您还可以看一看Numba，它很容易尝试，并且可能会稍微提高性能。

如您所见，这只是关于如何改进问题的一部分的提示。我敢打赌，还有其他几种方法可以进一步改进您的代码！

编辑

看起来Numba可能不是个坏主意。只需添加@numba.jit装饰器：

@numba.jit
def powernumba(x, n):
    result = [1., x]
    aux = x
    for i in range(2, n):
        aux *= x               
        result.append(aux)
    return result

然后：

In [52]: %timeit sum(power(1.1, 100))
100000 loops, best of 3: 7.67 µs per loop

In [51]: %timeit sum(powernumba(1.1, 100))
The slowest run took 5.64 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 2.64 µs per loop

看起来Numba可以在那里施展魔法。;-)

Answer 2

对于标量x

>>> import numpy as np
>>> x = 0.5
>>> jj = np.arange(10)
>>> x**jj
array([ 1.        ,  0.5       ,  0.25      ,  0.125     ,  0.0625    ,
        0.03125   ,  0.015625  ,  0.0078125 ,  0.00390625,  0.00195312])
>>> np.sin(jj*x)
array([ 0.        ,  0.47942554,  0.84147098,  0.99749499,  0.90929743,
        0.59847214,  0.14112001, -0.35078323, -0.7568025 , -0.97753012])
>>> (x**jj * np.sin(jj*x)).sum()
0.64489974041068521

注意numpy数组的sum方法的使用(等价地，使用np.sum而不是内置的sum)。

如果你的x本身是一个数组，使用广播：

>>> a = x[:, None]**jj
>>> a.shape
(3, 10)
>>> x[0]**jj == a[0]
array([ True,  True,  True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)

然后在第二个轴上求和：

>>> res = a * np.sin(jj * x[:, None])
>>> res.shape
(3, 10)
>>> res.sum(axis=1)
array([ 0.01230993,  0.0613201 ,  0.17154859])