7

我有一个向量在此处输入图像描述，并希望制作另一个相同长度的向量，其第 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 向量化更快？

4

3 回答 3

11

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

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

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

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

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

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

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

您可能可以使用 加速integrand（和减速）cython，但首先要专注于使用常规numpy代码获得最大的速度。

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

Withif w<0 w必须是标量。可以写成它与数组一起使用w吗？如果是这样，那么

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

可以改写为

 fn(x, ws).sum()

或者，由于x和w都是标量，因此使用math.expetc 而不是np.exp. log和相同abs。

我会尝试编写f(x,w)，因此它需要数组x和w，返回二维结果。如果是这样，那么tempandintegrand也可以与数组一起使用。由于quad提供了一个 scalar x，这在这里可能无济于事，但对于其他积分器来说，它可能会产生很大的不同。

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

于 2016-12-19T17:20:48.693 回答
0

您可以使用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]
于 2019-10-30T22:07:18.793 回答