问题简而言之
给定一个大的稀疏csr_matrix A和一个numpy数组B,构造numpy矩阵C的最快方法是什么,这样C [i,j] =所有k的sum(A [k,j]),其中B [k] ==我？
问题的细节
我找到了一个解决方案来做到这一点,但我并不满足于需要多长时间.我将首先解释问题,然后我的解决方案,然后显示我的代码,然后显示我的时间.
问题
我正在研究Python中的聚类算法,我想加快速度.我有一个稀疏的csr_matrix pam,其中每个人每篇文章他们购买了多少项目.此外,我有一个numpy数组聚类,其中表示该人所属的聚类.例：
pam pam.T clustering
article person
p [[1 0 0 0] a
e [0 2 0 0] r [[1 0 1 0 0 0] [0 0 0 0 1 1]
r [1 1 0 0] t [0 2 1 0 0 0]
s [0 0 1 0] i [0 0 0 1 0 1]
o [0 0 0 1] c [0 0 0 0 1 2]]
n [0 0 1 2]] l
e我想要计算的是acm：一个集群中所有人一起购买的物品数量.对于acm中的每列i,这相当于为pam.T添加了聚类[p] == i的p列.
acm
cluster
a
r [[2 0]
t [3 0]
i [1 1]
c [0 3]]
l
e解
首先,我创建另一个稀疏矩阵pcm,其中我指示每个元素[i,j],如果人i在群集j中.结果(当投射到密集矩阵时)：
pcm
cluster
p [[False True]
e [False True]
r [ True False]
s [False True]
o [False True]
n [ True False]]接下来,我用pcm将pam.T与pcm相乘以得到我想要的矩阵.
码
我编写了以下程序来测试这种方法在实践中的持续时间.
import numpy as np
from scipy.sparse.csr import csr_matrix
from timeit import timeit
def _clustering2pcm(clustering):
'''
Converts a clustering (np array) into a person-cluster matrix (pcm)
'''
N_persons = clustering.size
m_person = np.arange(N_persons)
clusters = np.unique(clustering)
N_clusters = clusters.size
m_data = [True] * N_persons
pcm = csr_matrix( (m_data, (m_person, clustering)), shape = (N_persons, N_clusters))
return pcm
def pam_clustering2acm():
'''
Convert a person-article matrix and a given clustering into an
article-cluster matrix
'''
global clustering
global pam
pcm = _clustering2pcm(clustering)
acm = csr_matrix.transpose(pam).dot(pcm).todense()
return acm
if __name__ == '__main__':
global clustering
global pam
N_persons = 200000
N_articles = 400
N_shoppings = 400000
N_clusters = 20
m_person = np.random.choice(np.arange(N_persons), size = N_shoppings, replace = True)
m_article = np.random.choice(np.arange(N_articles), size = N_shoppings, replace = True)
m_data = np.random.choice([1, 2], p = [0.99, 0.01], size = N_shoppings, replace = True)
pam = csr_matrix( (m_data, (m_person, m_article)), shape = (N_persons, N_articles))
clustering = np.random.choice(np.arange(N_clusters), size = N_persons, replace = True)
print timeit(pam_clustering2acm, number = 100)定时
事实证明,对于这100次运行,我需要5.1秒.其中3.6秒用于创建pcm.我觉得可以有一种更快的方法来计算这个矩阵而不创建一个临时的稀疏矩阵,但我没有看到没有循环的.有更快的建设方式吗？
编辑
在Martino的回答之后,我试图在集群和切片算法上实现循环,但这甚至更慢.现在需要12.5秒来计算acm 100次,如果我删除行acm [：,i] = pam [p,：].sum(axis = 0),则剩余4.1秒.
def pam_clustering2acm_loopoverclusters():
global clustering
global pam
N_articles = pam.shape[1]
clusters = np.unique(clustering)
N_clusters = clusters.size
acm = np.zeros([N_articles, N_clusters])
for i in clusters:
p = np.where(clustering == i)[0]
acm[:,i] = pam[p,:].sum(axis = 0)
return acm解决方法:
这比_clustering2pcm函数快50倍：
def pcm(clustering):
n = clustering.size
data = np.ones((n,), dtype=bool)
indptr = np.arange(n+1)
return csr_matrix((data, clustering, indptr))我没有查看源代码,但是当你传递CSR构造函数(data,(rows,cols))结构时,几乎可以肯定使用它来创建COO矩阵,然后将其转换为CSR.因为您的矩阵非常简单,所以很容易将实际的CSR矩阵描述数组放在一起,如上所述,并跳过所有这些.
这几乎可以将执行时间减少三个：
In [38]: %timeit pam_clustering2acm()
10 loops, best of 3: 36.9 ms per loop
In [40]: %timeit pam.T.dot(pcm(clustering)).A
100 loops, best of 3: 12.8 ms per loop
In [42]: np.all(pam.T.dot(pcm(clustering)).A == pam_clustering2acm())
Out[42]: True
