Python中2种稀疏矩阵的逐行乘法的特殊类型

Question

我正在寻找的:一种在Python中实现特殊乘法运算的方法,该矩阵恰好采用scipy稀疏(csr)格式.这是一种特殊的乘法,不是矩阵乘法,也不是Kronecker乘法,也不是Hadamard又是逐点乘法,并且在scipy.sparse中似乎没有任何内置支持.

所需操作:输出的每一行应包含两个输入矩阵中相应行元素的每个乘积的结果.所以开始具有两个相同尺寸的矩阵,每一个尺寸米由Ñ,结果应具有的尺寸米由N ^ 2.

它看起来像这样:

Python代码:

import scipy.sparse
A = scipy.sparse.csr_matrix(np.array([[1,2],[3,4]]))
B = scipy.sparse.csr_matrix(np.array([[0,5],[6,7]]))

# C is the desired product of A and B. It should look like:
C = scipy.sparse.csr_matrix(np.array([[0,5,0,10],[18,21,24,28]]))

什么是一个很好或有效的方法来做到这一点？我试过看看stackoverflow以及其他地方,到目前为止没有运气.到目前为止,听起来我最好的选择是在for循环中逐行操作,但这听起来很可怕,因为我的输入矩阵有几百万行和几千列,大多数是稀疏的.

Answer 1

在您的示例中,C是第一行和最后一行kron

In [4]: A=np.array([[1,2],[3,4]])
In [5]: B=np.array([[0,5],[6,7]])
In [6]: np.kron(A,B)
Out[6]: 
array([[ 0,  5,  0, 10],
       [ 6,  7, 12, 14],
       [ 0, 15,  0, 20],
       [18, 21, 24, 28]])
In [7]: np.kron(A,B)[[0,3],:]
Out[7]: 
array([[ 0,  5,  0, 10],
       [18, 21, 24, 28]])

kron包含相同的值np.outer,但它们的顺序不同.

对于大型密集阵列,einsum可能提供良好的速度:

np.einsum('ij,ik->ijk',A,B).reshape(A.shape[0],-1)

sparse.kron做同样的事情np.kron:

As = sparse.csr_matrix(A); Bs ...
sparse.kron(As,Bs).tocsr()[[0,3],:].A

sparse.kron 是用Python编写的,所以如果进行不必要的计算,你可能会修改它.

迭代解决方案似乎是:

sparse.vstack([sparse.kron(a,b) for a,b in zip(As,Bs)]).A

作为迭代,我不认为它比完全削减更快kron.但是,挖掘sparse.kron它的逻辑可能是我能做的最好的.

vstack使用bmat,所以计算是:

sparse.bmat([[sparse.kron(a,b)] for a,b in zip(As,Bs)])

但是bmat相当复杂,因此进一步简化这一点并不容易.

该np.einsum解决方案不能容易地扩展到稀疏-不存在sparse.einsum,以及所述中间产物是三维,其中稀疏不处理.

sparse.kron使用coo格式,这对于处理行没有好处.但是在这个函数的精神下工作,我已经找到了一个迭代csr格式矩阵行的函数.就像kron和bmat我建设data,row,col数组和构建coo_matrix从这些.这反过来可以转换为其他格式.

def test_iter(A, B):
    m,n1 = A.shape
    n2 = B.shape[1]
    Cshape = (m, n1*n2)
    data = np.empty((m,),dtype=object)
    col =  np.empty((m,),dtype=object)
    row =  np.empty((m,),dtype=object)
    for i,(a,b) in enumerate(zip(A, B)):
        data[i] = np.outer(a.data, b.data).flatten()
        #col1 = a.indices * np.arange(1,a.nnz+1) # wrong when a isn't dense
        col1 = a.indices * n2   # correction
        col[i] = (col1[:,None]+b.indices).flatten()
        row[i] = np.full((a.nnz*b.nnz,), i)
    data = np.concatenate(data)
    col = np.concatenate(col)
    row = np.concatenate(row)
    return sparse.coo_matrix((data,(row,col)),shape=Cshape)

使用这些小的2x2矩阵以及较大的矩阵(例如A1=sparse.rand(1000,2000).tocsr()),这比使用的版本快3倍bmat.对于足够大的矩阵,它比密集einsum版本(可能有内存错误)更好.