使用两个条件过滤Numpy数组

Question

我正在尝试运行自定义的kmeans聚类算法，但是在获取每个2-d numpy数组中每列（术语）的文档频率时遇到了问题。我的当前算法有两个numpy数组，一个原始数据集，按术语列出文档[2000L，9500L]，另一个是聚类分配[2000L，]。共有5个聚类。我需要做的是创建一个数组，列出每个聚类的文档频率 - 基本上是在每个列中计数，其中列号与不同数组中的行号匹配。输出将是一个[5L，9500L]数组（聚类x术语）。我无法找到一种方法来执行countif和group by的等效操作。以下是一些示例数据以及如果我仅使用2个聚类运行它时想要的输出：

import numpy as np

dataset = np.array[[1,2,0,3,0],[0,2,0,0,3],[4,5,2,3,0],[0,0,2,3,0]]
clusters = np.array[0,1,1,0]
#run code here to get documentFrequency
print documentFrequency
>> [1,1,1,2,0],[1,2,1,1,1]

我的想法是选择与每个聚类匹配的特定行，因为这样计数应该很容易。例如，如果我可以将数据拆分为以下数组：

cluster0 = np.array[[1,2,0,3,0],[0,0,2,3,0]]
cluster1 = np.array[[0,2,0,0,3],[4,5,2,3,0]]

任何指导或提示都将不胜感激！

- flyingmeatball

Answer 1

我认为没有简单的方法可以向量化你的代码，但是如果你只有几个聚类，你可以做一些显而易见的事情:

>>> cluster_count = np.max(clusters)+1
>>> doc_freq = np.zeros((cluster_count, dataset.shape[1]), dtype=dataset.dtype)
>>> for j in xrange(cluster_count):
...     doc_freq[j] = np.sum(dataset[clusters == j], axis=0)
... 
>>> doc_freq
array([[1, 2, 2, 6, 0],
       [4, 7, 2, 3, 3]])

- Jaime

---

谢谢你的建议，Jaime - 我会试一下。 - flyingmeatball

@flyingmeatball 这可能会对你有所帮助。我已经编写了 np.bincount 的 gufunc 版本，详见这里。如果你能够编译和安装它（如果你的系统已经正确设置，运行 python setup.py install 应该就可以了），那么你就可以像这样做：import new_gufuncs as ng; doc_freq = ng.bincount(clusters, dataset.T).T，所有循环都在 C 中进行。 - Jaime

Answer 2

正如@Jaime所说，如果你只有几个聚类，手动循环最小轴长度是有意义的常见技巧。通常这可以让你获得大部分矢量化的好处，而避免了聪明的头痛。

话虽如此，当你想要使用groupby时，你经常会发现自己处于一个需要高级工具的领域，比如pandas非常方便：

>>> pd.DataFrame(dataset).groupby(clusters).sum()
   0  1  2  3  4
0  1  2  2  6  0
1  4  7  2  3  3

如果需要的话，您可以轻松地回到一个 ndarray ：

>>> pd.DataFrame(dataset).groupby(clusters).sum().values
array([[1, 2, 2, 6, 0],
       [4, 7, 2, 3, 3]])

- DSM

---

谢谢，我会看一下Pandas - 你可能是对的，如果它不适合，最好尝试另一种方法。 - flyingmeatball

Answer 3

根据您的BLAS编译情况，将其编写为矩阵乘法可能会更快：

cvals = (clusters == np.arange(clusters.max()+1)[:,None]).astype(int)

cvals
array([[1, 0, 0, 1],
       [0, 1, 1, 0]])

np.dot(cvals,dataset)
array([[1, 2, 2, 6, 0],
       [4, 7, 2, 3, 3]])

让我们创建两个定义：

def loop(cvals,dataset):
     cluster_count = np.max(cvals)+1
     doc_freq = np.zeros((cluster_count, dataset.shape[1]), dtype=dataset.dtype)
     for j in xrange(cluster_count):
         doc_freq[j] = np.sum(dataset[cvals == j], axis=0)
     return doc_freq

def matrix_mult(clusters,dataset):
     cvals = (clusters == np.arange(clusters.max()+1)[:,None]).astype(dataset.dtype)
     return np.dot(cvals,dataset)

现在是一些时间：

arr = np.random.random((2000,9500))
cluster = np.random.randint(0,5,(2000))

np.allclose(loop(cluster,arr),matrix_mult(cluster,arr))
True

%timeit loop(cluster,arr)
1 loops, best of 3: 263 ms per loop

%timeit matrix_mult(cluster,arr)
100 loops, best of 3: 14.1 ms per loop

请注意，这是使用线程化的mkl BLAS。您的结果可能会有所不同。

- Daniel