各类特殊的Python

的2点稀疏矩阵的行由行乘法的，我找什么：Python中实现该恰巧在SciPy的稀疏（csr）格式矩阵特殊乘法运算的方式。这是一种特殊的乘法，而不是矩阵乘法也不Kronecker multiplication也不Hadamard又名pointwise乘，似乎没有任何内置在scipy.sparse支持。各类特殊的Python

所需的操作：所述输出的每一行中应包含对应的行中的两个输入矩阵的元素的每一个产品的结果。因此，与两个相同大小的矩阵，每一个尺寸米通过Ñ开始，结果应通过N^2具有尺寸米。

它看起来像这样：

graphic of desired multiplication operation

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循环中，但听起来可怕看到，因为我输入矩阵有几百万行和几千列，大多稀疏。

来源

2014-10-02 xenocyon

在你的榜样，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，而中间产品是3d，稀疏不处理。

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)

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

来源

2014-10-02 05:23:46 hpaulj

谢谢！我想到了运行kron并切分结果，但不确定这是否是一个好主意，因为在我的情况下，它会生成一个有数万亿行的矩阵。重写克朗是一个好主意。将介绍您建议的一些想法并查看最佳效果。 – xenocyon 2014-10-02 16:32:18

我已经完成了一个'kron'重写，与之前的'bmat'计算相比，这个重写提供了3倍的加速比。 – hpaulj 2014-10-03 00:14:20

谢谢！你重写的test_iter函数确实要快得多。我很好奇的一件事：如果输入矩阵中的行数增加，CPU时间增加比线性更差。这似乎很奇怪，因为行没有彼此的交互。任何想法，为什么这可能是这种情况？ – xenocyon 2014-10-03 20:41:08

非最佳的方式来做到这一点是每一行分别克朗：

def my_mult(A, B): 
    nrows = A.shape[0] 
    prodrows = [] 
    for i in xrange(0, nrows): 
     Arow = A.getrow(i) 
     Brow = B.getrow(i) 
     prodrow = scipy.sparse.kron(Arow,Brow) 
     prodrows.append(prodrow) 
    return scipy.sparse.vstack(prodrows)

这是大约3倍的性能比@ hpaulj的解决方案here更糟，因为可以通过运行下面的代码中可以看出：

A=scipy.sparse.rand(20000,1000, density=0.05).tocsr() 
B=scipy.sparse.rand(20000,1000, density=0.05).tocsr() 

# Check memory 
%memit C1 = test_iter(A,B) 
%memit C2 = my_mult(A,B) 

# Check time 
%timeit C1 = test_iter(A,B) 
%timeit C2 = my_mult(A,B) 

# Last but not least, check correctness! 
print (C1 - C2).nnz == 0

结果：

hpaulj的方法：

peak memory: 1993.93 MiB, increment: 1883.80 MiB 
1 loops, best of 3: 6.42 s per loop

这种方法：

peak memory: 2456.75 MiB, increment: 1558.78 MiB 
1 loops, best of 3: 18.9 s per loop

来源

2014-10-15 01:46:20 xenocyon

各类特殊的Python

回答

相关问题