Numpy数组过滤两个标准答案

【问题标题】：Numpy array filtering by two criteriaNumpy数组过滤两个标准
【发布时间】：2013-11-19 16:58:03
【问题描述】：

我正在尝试运行自定义 kmeans 聚类算法，但无法按聚类获取二维 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]]

任何方向或指针将不胜感激！

【问题讨论】：

标签： python numpy cluster-analysis k-means

【解决方案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 的建议 - 我会尝试一下。
@flyingmeatball 您可能对此感兴趣。我写了一个 gufunc 版本的 np.bincount，见 here。如果您可以编译和安装它（如果您的系统设置正确，则运行 python setup.py install 应该可以解决问题），然后您可以执行类似 import new_gufuncs as ng; doc_freq = ng.bincount(clusters, dataset.T).T 的操作，所有循环都在 C 中进行。

【解决方案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]])

【讨论】：

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

【解决方案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。您的里程会有所不同。

【讨论】：