我有一个稀疏的csc矩阵,其中许多元素为零,我想要计算每行所有列元素的乘积。
即:
A = [[1,2,0,0],
[2,0,3,0]]
应该被转换为:
V = [[2,
6]]
使用numpy的密集矩阵,可以通过将所有零值替换为一值并使用
A.prod(1)
来实现。但是,由于密集矩阵太大,这不是一个选择。是否有任何方法可以在不将稀疏矩阵转换为密集矩阵的情况下完成此操作?我有一个稀疏的csc矩阵,其中许多元素为零,我想要计算每行所有列元素的乘积。
即:
A = [[1,2,0,0],
[2,0,3,0]]
应该被转换为:
V = [[2,
6]]
A.prod(1)
来实现。但是,由于密集矩阵太大,这不是一个选择。是否有任何方法可以在不将稀疏矩阵转换为密集矩阵的情况下完成此操作?np.multiply.reduceat
对这些元素的相应值进行乘法运算,以获得所需的输出。from scipy import sparse
from scipy.sparse import csc_matrix
r,c,v = sparse.find(a) # a is input sparse matrix
out = np.zeros(a.shape[0],dtype=a.dtype)
unqr, shift_idx = np.unique(r,return_index=1)
out[unqr] = np.multiply.reduceat(v, shift_idx)
样例运行 -
In [89]: # Let's create a sample csc_matrix
...: A = np.array([[-1,2,0,0],[0,0,0,0],[2,0,3,0],[4,5,6,0],[1,9,0,2]])
...: a = csc_matrix(A)
...:
In [90]: a
Out[90]:
<5x4 sparse matrix of type '<type 'numpy.int64'>'
with 10 stored elements in Compressed Sparse Column format>
In [91]: a.toarray()
Out[91]:
array([[-1, 2, 0, 0],
[ 0, 0, 0, 0],
[ 2, 0, 3, 0],
[ 4, 5, 6, 0],
[ 1, 9, 0, 2]])
In [92]: out
Out[92]: array([ -2, 0, 6, 120, 0, 18])
方法 #2: 我们正在执行基于二进制的乘法。我们有一个基于二进制的求和解决方案,使用 np.bincount
。因此,一个可以使用的技巧是将数字转换为对数数字,进行基于二进制的求和,然后使用 指数函数
(log的反函数)将其转换回原始格式,就这样!对于负数,我们可能需要添加一步或更多步骤,但让我们看看非负数的实现情况 -
r,c,v = sparse.find(a)
out = np.exp(np.bincount(r,np.log(v),minlength = a.shape[0]))
out[np.setdiff1d(np.arange(a.shape[0]),r)] = 0
使用非负数进行的样本运行 -
In [118]: a.toarray()
Out[118]:
array([[1, 2, 0, 0],
[0, 0, 0, 0],
[2, 0, 3, 0],
[4, 5, 6, 0],
[1, 9, 0, 2]])
In [120]: out # Using listed code
Out[120]: array([ 2., 0., 6., 120., 18.])
In [51]: A=np.array([[1,2,0,0],[0,0,0,0],[2,0,3,0]])
In [52]: M=sparse.csr_matrix(A)
In [56]: Ml=M.tolil()
In [57]: Ml.data
Out[57]: array([[1, 2], [], [2, 3]], dtype=object)
取每个数的乘积:
In [58]: np.array([np.prod(i) for i in Ml.data])
Out[58]: array([ 2., 1., 6.])
csr
格式中,值存储为:In [53]: M.data
Out[53]: array([1, 2, 2, 3], dtype=int32)
In [54]: M.indices
Out[54]: array([0, 1, 0, 2], dtype=int32)
In [55]: M.indptr
Out[55]: array([0, 2, 2, 4], dtype=int32)
indptr
给出行值的起始位置。在 csr
(和 csc
) 矩阵上进行计算的代码通常会执行类似这样的计算(但已编译):
In [94]: lst=[]; i=M.indptr[0]
In [95]: for j in M.indptr[1:]:
...: lst.append(np.product(M.data[i:j]))
...: i = j
In [96]: lst
Out[96]: [2, 1, 6]
In [137]: M.A
Out[137]:
array([[-1, 2, 0, 0],
[ 0, 0, 0, 0],
[ 2, 0, 3, 0],
[ 4, 5, 6, 0],
[ 1, 9, 0, 2]], dtype=int32)
In [138]: foo(M)
Out[138]: [-2, 1, 6, 120, 18]
unique
和reduceat
In [139]: divk(M)
Out[139]: array([ -2, 0, 6, 120, 18], dtype=int32)
indptr
的Reduceat很简单:In [140]: np.multiply.reduceat(M.data,M.indptr[:-1])
Out[140]: array([ -2, 2, 6, 120, 18], dtype=int32)
indptr
值为[2,2,...],reduceat
将使用M.data[2]
)。def wptr(M, empty_val=1):
res = np.multiply.reduceat(M.data, M.indptr[:-1])
mask = np.diff(M.indptr)==0
res[mask] = empty_val
return res
使用更大的矩阵
Mb=sparse.random(1000,1000,.1,format='csr')
wptr
比 Divaker 版本快了大约 30 倍。import numpy as np
print [[np.prod([x for x in A[i] if x!=0 ]) for i in range(len(A))]]
reduceat
中使用csr
的indptr
会更快。 - hpauljfloat128
值,而bincount
还不支持(但可能会在未来支持)。我在使用你的第一种方法时遇到了一些问题,直到我意识到,在我的安装中,find
的结果不是按行排序的,而是按列排序的(我像你一样使用了csc_matrix
)。在弄清楚之后,我能够使用转置矩阵解决它。 - Martin