9

我有一个 numpy 大小的数组

arr.size = (200, 600, 20).

我想计算scipy.stats.kendalltau最后两个维度的每个成对组合。例如：

kendalltau(arr[:, 0, 0], arr[:, 1, 0])
kendalltau(arr[:, 0, 0], arr[:, 1, 1])
kendalltau(arr[:, 0, 0], arr[:, 1, 2])
...
kendalltau(arr[:, 0, 0], arr[:, 2, 0])
kendalltau(arr[:, 0, 0], arr[:, 2, 1])
kendalltau(arr[:, 0, 0], arr[:, 2, 2])
...
...
kendalltau(arr[:, 598, 20], arr[:, 599, 20])

这样我就涵盖arr[:, i, xi]了 witharr[:, j, xj]和i < j,xi in [0,20)的所有组合xj in [0, 20)。这是(600 choose 2) * 400单独的计算，但由于每个计算都需要0.002 s在我的机器上进行，因此使用多处理模块所花费的时间不会超过一天。

遍历这些列（使用i<j）的最佳方法是什么？我想我应该避免类似的事情

for i in range(600):
    for j in range(i+1, 600):
        for xi in range(20):
            for xj in range(20):

这样做的最简单的方法是什么？

编辑：我更改了标题，因为 Kendall Tau 对这个问题并不重要。我意识到我也可以做类似的事情

import itertools as it
for i, j in it.combinations(xrange(600), 2):
    for xi, xj in product(xrange(20), xrange(20)):

但是 numpy 必须有一种更好、更矢量化的方式。

4

2 回答 2

14

将此类向量化的一般方法是使用广播来创建集合与自身的笛卡尔积。arr在您的情况下，您有一个shape数组(200, 600, 20)，因此您可以对其进行两种查看：

arr_x = arr[:, :, np.newaxis, np.newaxis, :] # shape (200, 600, 1, 1, 20)
arr_y = arr[np.newaxis, np.newaxis, :, :, :] # shape (1, 1, 200, 600, 20)

为清楚起见，以上两行已被扩展，但我通常会写成等价的：

arr_x = arr[:, :, None, None]
arr_y = arr

如果您有一个矢量化函数 ，f它在除最后一个维度之外的所有维度上进行广播，那么您可以执行以下操作：

out = f(arr[:, :, None, None], arr)

然后out将是一个 shape 数组(200, 600, 200, 600)，其中包含out[i, j, k, l]的值f(arr[i, j], arr[k, l])。例如，如果你想计算所有的成对内积，你可以这样做：

from numpy.core.umath_tests import inner1d

out = inner1d(arr[:, :, None, None], arr)

不幸scipy.stats.kendalltau的是没有像这样矢量化。根据文档

“如果数组不是一维的，它们将被展平为一维。”

所以你不能像这样去做，你最终会做 Python 嵌套循环，无论是明确地写出来，使用itertools还是伪装它np.vectorize。这会很慢，因为 Python 变量的迭代，并且因为每个迭代步骤都有一个 Python 函数，这都是昂贵的操作。

请注意，当您可以采用矢量化方式时，有一个明显的缺点：如果您的函数是可交换的，即 if f(a, b) == f(b, a)，那么您需要进行两倍的计算。根据您的实际计算成本，这通常会被没有任何 Python 循环或函数调用所带来的速度增加所抵消。

于 2013-08-09T21:18:08.357 回答
0

如果您不想使用递归，通常应该使用itertools.combinations。 没有具体原因（afaik）为什么这会导致您的代码运行速度变慢。计算密集型部分仍由 numpy 处理。Itertools 还具有可读性的优势。

于 2013-08-09T20:39:21.030 回答