NumPy矢量化与集成

Question

我有一个矢量 并希望制作另一个长度相同的矢量,其第k个分量是

问题是:我们如何才能将其矢量化以提高速度？NumPy vectorize()实际上是for循环,所以它不计算.

Veedrac指出" 没有办法将纯Python函数应用于NumPy数组的每个元素而不会多次调用它 ".由于我使用的是NumPy函数而不是"纯Python"函数,我认为它可以进行矢量化,但我不知道如何.

import numpy as np
from scipy.integrate import quad
ws = 2 * np.random.random(10) - 1
n  = len(ws)
integrals = np.empty(n)

def f(x, w):
    if w < 0: return np.abs(x * w)
    else:     return np.exp(x) * w

def temp(x): return np.array([f(x, w) for w in ws]).sum()

def integrand(x, w): return f(x, w) * np.log(temp(x))

## Python for loop
for k in range(n):
    integrals[k] = quad(integrand, -1, 1, args = ws[k])[0]

## NumPy vectorize
integrals = np.vectorize(quad)(integrand, -1, 1, args = ws)[0]

另外,Cython for循环总是比NumPy矢量化更快吗？

Answer 1

该函数quad执行自适应算法,这意味着它执行的计算取决于所集成的特定事物.这原则上不能矢量化.

在您的情况下,for长度为10 的循环不是问题.如果程序需要很长时间,那是因为集成需要很长时间,而不是因为你有一个for循环.

当您绝对需要向量化矢量化时(不在上面的示例中),请使用非自适应方法,并了解精度可能会受到影响.这些可以直接应用于通过评估某些规则间隔的1D阵列(a linspace)上的所有函数而获得的2D NumPy数组.您必须自己选择linspace,因为这些方法不是自适应的.

• numpy.trapz是最简单和最不精确的
• scipy.integrate.simps同样易于使用且更精确(Simpson的规则需要奇数个样本,但该方法也可以使用偶数).
• scipy.integrate.romb原则上比Simpson具有更高的精度(对于平滑数据),但它要求样本数量2**n+1为某个整数n.