如何将 for 循环下推到 numpy答案

【问题标题】：How to push the for-loop down to numpy如何将 for 循环下推到 numpy
【发布时间】：2013-09-28 06:40:53
【问题描述】：

下面的代码完全符合我的要求（它是克里金法的一部分）。但问题是它运行得太慢了，我想知道是否有任何选项可以将 for 循环推到 numpy？如果我推出 numpy.sum，并在那里使用轴参数，它会加快一点，但显然这不是瓶颈。关于如何将 forloop 下推到 numpy 以加快速度或其他方式来加快速度的任何想法？）

# n = 2116
print GRZVV.shape  # (16309, 2116)
print GinvVV.shape  # (2117, 2117) 
VVg = numpy.empty((GRZVV.shape[0]))

for k in xrange(GRZVV.shape[0]):
    GRVV = numpy.empty((n+1, 1))
    GRVV[n, 0] = 1
    GRVV[:n, 0] = GRZVV[k, :]
    EVV = numpy.array(GinvVV * GRVV)  # GinvVV is numpy.matrix
    VVg[k] = numpy.sum(EVV[:n, 0] * VV)

我发布了 ndarrays n 矩阵的维度来清除一些东西

编辑：VV 的形状是 2116

【问题讨论】：

如果VV.shape == (16309,)，你怎么能把它乘以EVV[:n, 0]，它的形状是(n,)？
也许循环的最后一行应该有EVV[:n, 0] * VV[k]，这似乎是@Jaime 的回答所假设的。
@askewchan 的形状是 2116，那里有点混乱，我编辑了它，

标签： python numpy kriging

【解决方案1】：

您可以执行以下操作来代替 k 上的循环（运行时 ~3 秒）：

tmp = np.concatenate((GRZVV, np.ones((16309,1),dtype=np.double)), axis=1)
EVV1 = np.dot(GinvVV, tmp.T)
#Changed line below based on *askewchan's* recommendation
VVg1 = np.sum(np.multiply(EVV1[:n,:],VV[:,np.newaxis]), axis=0)

【讨论】：

这给出了与@usethedeathstar 的代码相同的结果，并且在我的机器上运行速度提高了 15 倍。
不需要磁贴调用，因为np.multiply 广播。将其更改为：VVg1 = np.sum(np.multiply(EVV1[:n,:],VV[:,np.newaxis]), axis=0) 以加快速度。
+1 对于大型数组，比 np.einsum 快得多，不知道为什么……
小问题：VVg1 仍然是 numpy.matrix 而不是 ndarray？这会影响 np.multiply 吗？我假设/希望不是，或者我应该在 np.dot 之后对 EVV1 上的 numpy.array 进行强制转换？
另一个问题： [:,numpy.newaxis] 和使用 [:,None] 有什么区别？为什么 numpy.multiply 而不是只使用 "*" ？并且加速是 12 倍（这是 NICE）

【解决方案2】：

您基本上是取GRZVV 的每一行，在末尾附加一个 1，将其与GinvVV 相乘，然后将向量中的所有元素相加。如果你没有做“追加 1”的事情，你可以在没有循环的情况下完成这一切：

VVg = np.sum(np.dot(GinvVV[:, :-1], GRZVV.T), axis=-1) * VV

甚至：

VVg = np.einsum('ij,kj->k', GinvVV[:, :-1], GRZVV) * VV

我们如何处理额外的 1？好吧，来自矩阵乘法的结果向量将增加GinvVV[:, -1] 中的相应值，当您将它们全部相加时，该值将增加np.sum(GinvVV[:, -1])。所以我们可以简单地计算一次并将其添加到返回向量中的所有项目中：

VVg = (np.einsum('ij,kj->k', GinvVV[:-1, :-1], GRZVV) + np.sum(GinvVV[:-1, -1])) * VV

如果VV 是标量，则上述代码有效。如果它是一个形状为(n,) 的数组，那么以下将起作用：

GinvVV = np.asarray(GinvVV)
VVgbis = (np.einsum('ij,kj->k', GinvVV[:-1, :-1]*VV[:, None], GRZVV) +
          np.dot(GinvVV[:-1, -1], VV))

【讨论】：

在 @usethedeathstar 的编辑下，VV.shape 现在是 2116，因此不会在您的解决方案中广播（因为 Wg.shape 是 16309）
认为它是一个标量。使用完整的数组，直到最后才折叠数组，如 Joel 的回答，比上面的大型数组快得多。