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

所需操作：输出的每一行应包含两个输入矩阵中相应行元素的每个乘积的结果.因此,从两个大小相等的矩阵开始,每个矩阵的维数为n,结果的维度应为n乘以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循环中逐行操作,但这听起来很可怕,因为我的输入矩阵有几百万行和几千列,大多数是稀疏的.

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

作为迭代,我不认为它比削减完整的克朗更快.但是,如果没有深入挖掘sparse.kron的逻辑,它可能是我能做的最好的.

vstack使用bmat,因此计算如下：

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

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

np.einsum解决方案不能轻易扩展为稀疏 – 没有sparse.einsum,中间产品是3d,稀疏不能处理.

sparse.kron使用了coo格式,这对于处理行是没有用的.但是按照这个函数的精神工作,我已经找到了一个迭代csr格式矩阵行的函数.就像kron和bmat一样,我正在构造数据,行,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)

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