NumPy矢量化与集成 [英] NumPy vectorization with integration

问题描述

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

I have a vector and wish to make another vector of the same length whose k-th component is

问题是:我们如何向量化它以提高速度? NumPy vectorize()实际上是一个for循环，因此不计算在内.

The question is: how can we vectorize this for speed? NumPy vectorize() is actually a for loop, so it doesn't count.

Veedrac指出"如果不多次调用NumPy数组，就无法将纯Python函数应用于NumPy数组的每个元素".由于我使用的是NumPy函数，而不是纯Python"函数，因此我认为可以进行矢量化，但是我不知道该怎么做.

Veedrac pointed out that "There is no way to apply a pure Python function to every element of a NumPy array without calling it that many times". Since I'm using NumPy functions rather than "pure Python" ones, I suppose it's possible to vectorize, but I don't know how.

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 loop是否总是比NumPy向量化更快?

On a side note, is a Cython for loop always faster than NumPy vectorization?

推荐答案

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

The function quad executes an adaptive algorithm, which means the computations it performs depend on the specific thing being integrated. This cannot be vectorized in principle.

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

In your case, a for loop of length 10 is a non-issue. If the program takes long, it's because integration takes long, not because you have a for loop.

当您绝对需要向量化积分时(在上面的示例中不是)，请使用非自适应方法，并要注意精度可能会受到影响.这些可以直接应用于2D NumPy数组，该数组是通过在规则排列的1D数组(a linspace)上评估所有函数而获得的.由于方法不具有自适应性，因此您必须自己选择linspace.

When you absolutely need to vectorize integration (not in the example above), use a non-adaptive method, with the understanding that precision may suffer. These can be directly applied to a 2D NumPy array obtained by evaluating all of your functions on some regularly spaced 1D array (a linspace). You'll have to choose the linspace yourself since the methods aren't adaptive.

• numpy.trapz 是最简单，最少的精确
• scipy.integrate.simps 是同样易于使用且更精确(辛普森法则要求样本数量为奇数，但该方法也可以使用偶数).
• scipy.integrate.原则上，romb 比Simpson更高(对于平滑数据)，但是对于某些整数n，它要求样本数为2**n+1.

• numpy.trapz is the simplest and least precise
• scipy.integrate.simps is equally easy to use and more precise (Simpson's rule requires an odd number of samples, but the method works around having an even number, too).
• scipy.integrate.romb is in principle of higher accuracy than Simpson (for smooth data) but it requires the number of samples to be 2**n+1 for some integer n.

这篇关于NumPy矢量化与集成的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！