带有集成的 NumPy 矢量化

Question

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

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

维德拉克指出“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 执行自适应算法，这意味着它执行的计算取决于所集成的特定事物。这原则上不能向量化。

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

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

• numpy.trapz 是最简单最不精确的
• scipy.integrate.simps 同样易于使用且更精确（辛普森规则需要奇数个样本，但该方法也适用于偶数个样本）。
• scipy.integrate.romb 原则上比 Simpson 具有更高的准确度（对于平滑数据），但对于某个整数 n，它要求样本数为 2**n+1。

Answer 2

@zaq's 关注quad 的答案是正确的。所以我会看看这个问题的其他一些方面。

在最近的https://stackoverflow.com/a/41205930/901925 中，我认为vectorize 在您需要将完整的广播机制应用于只接受标量值的函数时最有价值。您的 quad 有资格接受标量输入。但是您只在一个数组ws 上进行迭代。传递给您的函数的x 是由quad 本身生成的。 quad 和 integrand 仍然是 Python 函数，即使它们使用 numpy 操作。

cython 改进了低级迭代，它可以转换为C 代码。您的主要迭代处于较高级别，调用导入的函数quad。 Cython 无法触摸或重写它。

您也许可以使用 cython 加速 integrand（并继续减速），但首先要专注于使用常规 numpy 代码获得最大速度。

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

if w<0 w 必须是标量。它可以写成它与数组w一起工作吗？如果是的话，那么

 np.array([f(x, w) for w in ws]).sum()

可以改写为

 fn(x, ws).sum()

另外，由于x 和w 都是标量，您可以通过使用math.exp 等而不是np.exp 来提高速度。 log 和 abs 相同。

我会尝试编写f(x,w)，因此它需要x 和w 的数组，返回二维结果。如果是这样，那么temp 和integrand 也适用于数组。由于quad 提供了一个标量x，这在这里可能无济于事，但对于其他集成商来说，它可能会产生很大的不同。

如果 f(x,w) 可以在 x=np.linspace(-1,1,n) 和 ws 的常规 nx10 网格上进行评估，那么（某种类型的）积分只需要在该空间上进行几次求和。

Answer 3

您可以使用quadpy 进行完全矢量化计算。您必须首先调整您的函数以允许矢量输入，但这很容易完成：

import numpy as np
import quadpy

np.random.seed(0)
ws = 2 * np.random.random(10) - 1


def f(x):
    out = np.empty((len(ws), *x.shape))
    out0 = np.abs(np.multiply.outer(ws, x))
    out1 = np.multiply.outer(ws, np.exp(x))
    out[ws < 0] = out0[ws < 0]
    out[ws >= 0] = out1[ws >= 0]
    return out


def integrand(x):
    return f(x) * np.log(np.sum(f(x), axis=0))


val, err = quadpy.quad(integrand, -1, +1, epsabs=1.0e-10)
print(val)

[0.3266534  1.44001826 0.68767868 0.30035222 0.18011948 0.97630376
 0.14724906 2.62169217 3.10276876 0.27499376]