稀疏矩阵列元素相乘

3

我有一个稀疏的csc矩阵,其中许多元素为零,我想要计算每行所有列元素的乘积。

即:

 A = [[1,2,0,0],
      [2,0,3,0]]

应该被转换为:

V = [[2,
      6]]

使用numpy的密集矩阵,可以通过将所有零值替换为一值并使用A.prod(1)来实现。但是,由于密集矩阵太大,这不是一个选择。是否有任何方法可以在不将稀疏矩阵转换为密集矩阵的情况下完成此操作?
3个回答

2
方法一:我们可以使用稀疏元素的行索引作为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: 我们正在执行基于二进制的乘法。我们有一个基于二进制的求和解决方案,使用 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.])

reduceat 中使用 csrindptr 会更快。 - hpaulj
我查看了你的两种方法,发现它们非常有帮助。我开始使用对数方法,这个方法很好用,直到我不得不放弃它,因为我需要使用 float128 值,而 bincount 还不支持(但可能会在未来支持)。我在使用你的第一种方法时遇到了一些问题,直到我意识到,在我的安装中,find 的结果不是按行排序的,而是按列排序的(我像你一样使用了 csc_matrix)。在弄清楚之后,我能够使用转置矩阵解决它。 - Martin

1
制作一个样本:
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]
Divakar的代码,使用了uniquereduceat
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 倍。
关于在稀疏矩阵的行之间计算值的更多讨论: Scipy.sparse.csr_matrix: 如何获取前十个值和索引?

我对稀疏矩阵不是很了解,但这看起来很有前途! - Divakar

0
你可以使用 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))]]

请使用[编辑]链接来解释这段代码的工作原理,而不仅仅是提供代码,因为解释更有可能帮助未来的读者。 - Jed Fox

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接