带集成的 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 函数而不是 "pure 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 循环。

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

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

Answer 2

@zaq's 专注于 quad 的回答是正确的。所以我会看看问题的其他方面。

最近我认为当您需要将完整的广播机制应用于仅采用标量值的函数时，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

With 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 网格上计算，那么一个积分（某种）只需要对 space.

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]

带集成的 NumPy 矢量化

NumPy vectorization with integration

numpy

vectorization

quad