python - 如何以矢量化形式编写此 numpy 代码？

Question

我在 python 中有以下函数，我无法弄清楚如何以矢量化形式表达。对我来说，cov 是一个形状为 (2,2) 的 numpy 数组，mu 是形状为 (2,) 的平均向量，xtp 是形状为 (~50000,2)。我知道 scipy 提供 scipy.stats.multivariate_normal.pdf 但我正在尝试学习如何编写高效的矢量化代码。请

def mvnpdf(xtp, mu, cov):
    temp = np.zeros(xtp.shape[0])
    i = 0
    length = xtp.shape[0]
    const = 1 / ( ((2* np.pi)**(len(mu)/2)) * (np.linalg.det(cov)**(1/2)) )
    inv = np.linalg.inv(cov)
    while i < length:
        x = xtp[i]-mu
        exponent = (-1/2) * (x.dot(inv).dot(x))
        temp[i] =  (const * np.exp(exponent))
        i+=1
    return temp

Answer 1

矢量化唯一棘手的部分是双 .dot。让我们隔离一下:

x = xtp - mu  # move this out of the loop
ddot = [i.dot(inv).dot(i) for i in x]
temp = const * np.exp(-0.5 * ddot)

将其放入您的代码中，看看它是否产生相同的结果。

“矢量化”点有几种方法。我首先喜欢尝试的是 einsum。在我的测试中，这是等效的:

ddot = np.einsum('ij,jk,ik->i',x,inv,x)

我建议尝试一下，看看它是否有效并加快速度。并在交互式 shell 中使用较小的数组(不是 ~50000)进行这些计算。

我正在测试

In [225]: x
Out[225]: 
array([[  0.,   2.],
       [  1.,   3.],
       ...
       [  7.,   9.],
       [  8.,  10.],
       [  9.,  11.]])
In [226]: inv
Out[226]: 
array([[ 1.,  0.],
       [ 0.,  1.]])

由于这是一个学习练习，我会将详细信息留给您。

使用 (2,2)，计算一个 cov 可能会更快，因为它是明确完成的，而不是使用 det 和 inv 函数。但消耗时间的正是 length 迭代。