带集成的 NumPy 矢量化

Question

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



问题是:我们如何将其矢量化以提高速度？ NumPy vectorize() 实际上是一个 for 循环，所以它不算数。

Veedrac 指出“There is no way to apply a pure Python function to every element of a NumPy array without calling it that many times”。由于我使用的是 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环形。

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

numpy.trapz是最简单和最不精确的

scipy.integrate.simps同样易于使用且更精确(辛普森规则需要奇数个样本，但该方法也适用于偶数个样本)。

scipy.integrate.romb原则上比Simpson(对于平滑数据)精度更高，但它要求样本数为2**n+1对于一些整数 n .