稀疏矩阵的乘积列元素

Question

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

即：

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

应转换为：

V = [[2, 
     6]] 

使用numpy的密集矩阵，这可以通过与一个值替换所有零个值，并使用A.prod(1)来完成。然而，这不是一种选择，因为密集的矩阵太大了。

是否有任何方法可以在不将稀疏矩阵转换为密集矩阵的情况下完成此任务？

Answer 1

方法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：我们执行的是基于斌乘法。我们有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.]) 

Answer 2

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

Answer 3

做一个样本：

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的代码unique和reduceat

In [139]: divk(M) 
Out[139]: array([ -2, 0, 6, 120, 18], dtype=int32) 

（不同的空行值）。

Reduceat与indptr仅仅是：需要

In [140]: np.multiply.reduceat(M.data,M.indptr[:-1]) 
Out[140]: array([ -2, 2, 6, 120, 18], dtype=int32) 

为空第二行中的值是固定的（与[2,2，...]的indptr值，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大约是30倍比Divaker的版本快。

横跨稀疏矩阵的行计算值的更多讨论： Scipy.sparse.csr_matrix: How to get top ten values and indices?