我有一个稀疏的 csc 矩阵，其中包含许多零元素，我想计算每行所有列元素的乘积。
即：
A = [[1,2,0,0],
[2,0,3,0]]
应转换为：
V = [[2,
6]]
使用 numpy 密集矩阵可以通过将所有零值替换为一个值并使用 A.prod(1) 来完成。然而，这不是一个选择，因为密集矩阵会太大。
有什么方法可以在不将稀疏矩阵转换为密集矩阵的情况下实现这一点？
【问题讨论】：
标签： python numpy matrix scipy sparse-matrix
方法#1：我们可以使用稀疏元素的行索引作为 ID，并将这些元素的相应值与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： 我们正在执行基于 bin 的乘法。我们有np.bincount 的基于bin 的求和解决方案。因此，可以在这里使用的一个技巧是将数字转换为对数，执行基于 bin 的求和，然后使用exponential（对数的反向）转换回原始格式，就是这样！对于负数，我们可能会添加一个或更多的步骤，但让我们看看非负数的实现是什么样的 -
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.])
【讨论】：
reduceat 中使用csr indptr 甚至更快。
bincount 不支持的 float128 值（还）。我对您的第一种方法遇到了一些麻烦，直到我意识到至少在我的安装中find 的结果不是按行而是按列排序（我像您一样使用了csc_matrix）。在弄清楚之后，我能够使用转置矩阵来解决它。
您可以使用 numpy 模块中的 prod() 方法计算 A 的每个子列表中所有元素的乘积，同时不考虑值为 0 的元素。
import numpy as np
print [[np.prod([x for x in A[i] if x!=0 ]) for i in range(len(A))]]
【讨论】：
制作样本：
In [51]: A=np.array([[1,2,0,0],[0,0,0,0],[2,0,3,0]])
In [52]: M=sparse.csr_matrix(A)
在lil 格式中，每一行的值都存储在一个列表中。
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]
使用Diavaker的测试矩阵：
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 的Divakar 代码
In [139]: divk(M)
Out[139]: array([ -2, 0, 6, 120, 18], dtype=int32)
（空行的不同值）。
使用indptr 减少操作很简单：
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 倍。
更多关于计算稀疏矩阵行的值的讨论： Scipy.sparse.csr_matrix: How to get top ten values and indices?
【讨论】：