将卷积应用于scipy.sparse矩阵

如何解决将卷积应用于scipy.sparse矩阵

我尝试在scipy.sparse矩阵上计算卷积。这是代码：

import numpy as np
import scipy.sparse,scipy.signal

M = scipy.sparse.csr_matrix([[0,1,0],[1,1],[0,0]])
kernel = np.ones((3,3))
kernel[1,1]=0
X = scipy.signal.convolve(M,kernel,mode='same')

哪个会产生以下错误：

ValueError: volume and kernel should have the same dimensionality

计算scipy.signal.convolve(M.todense(),mode='same')提供了预期的结果。但是，我想保持计算稀疏。

更笼统地说，我的目标是计算稀疏矩阵M的1跳邻域和。如果您有什么好方法，如何在稀疏矩阵上进行计算，我很乐意听到！

编辑：

我刚刚尝试了一种针对特定内核（邻居之和）的解决方案，该解决方案的速度实际上并不比密集版本快（不过，我并未尝试过高维度）。这是代码：

row_ind,col_ind = M.nonzero() 
X = scipy.sparse.csr_matrix((M.shape[0]+2,M.shape[1]+2))
for i in [0,2]:
    for j in [0,2]:
        if i!= 1 or j !=1:
            X += scipy.sparse.csr_matrix( (M.data,(row_ind+i,col_ind+j)),(M.shape[0]+2,M.shape[1]+2))
X = X[1:-1,1:-1]

解决方法

In [1]: from scipy import sparse,signal
In [2]: M = sparse.csr_matrix([[0,1,0],[1,1],[0,0]])
   ...: kernel = np.ones((3,3))
   ...: kernel[1,1]=0
In [3]: X = signal.convolve(M.A,kernel,mode='same')
In [4]: X
Out[4]: 
array([[2.,1.,2.,1.],[2.,4.,3.,[1.,2.],1.]])

为什么海报显示可运行的代码，而不显示结果？我们大多数人无法像这样运行代码。

In [5]: M.A
Out[5]: 
array([[0,0]])

您的替代方法-结果为稀疏矩阵时，将填充所有值。即使M较大且较稀疏，X也会较密集。

In [7]: row_ind,col_ind = M.nonzero()
   ...: X = sparse.csr_matrix((M.shape[0]+2,M.shape[1]+2))
   ...: for i in [0,2]:
   ...:     for j in [0,2]:
   ...:         if i!= 1 or j !=1:
   ...:             X += sparse.csr_matrix( (M.data,(row_ind+i,col_ind+j)),(M
   ...: .shape[0]+2,M.shape[1]+2))
   ...: X = X[1:-1,1:-1]
In [8]: X
Out[8]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 16 stored elements in Compressed Sparse Row format>
In [9]: X.A
Out[9]: 
array([[2.,1.]])

这里是构建coo样式输入并仅在末尾形成矩阵的替代方法。请记住，将重复的坐标相加。这在FEM刚度矩阵构造中很方便，并且在这里也很合适。

In [10]: row_ind,col_ind = M.nonzero()
    ...: data,row,col = [],[],[]
    ...: for i in [0,2]:
    ...:     for j in [0,2]:
    ...:         if i!= 1 or j !=1:
    ...:             data.extend(M.data)
    ...:             row.extend(row_ind+i)
    ...:             col.extend(col_ind+j)
    ...: X = sparse.csr_matrix( (data,(row,col)),(M.shape[0]+2,M.shape[1]+2)
    ...: )
    ...: X = X[1:-1,1:-1]
In [11]: X
Out[11]: 
<4x4 sparse matrix of type '<class 'numpy.int64'>'
    with 16 stored elements in Compressed Sparse Row format>
In [12]: X.A
Out[12]: 
array([[2,2,[2,4,3,2],1]])

===

我的方法明显更快（但仍远远落后于密集卷积）。 sparse.csr_matrix(...)的运行速度很慢，因此重复执行并不是一个好主意。而且稀疏添加也不是很好。

In [13]: %%timeit
    ...: row_ind,1:-1]
    ...: 
    ...: 
793 µs ± 20 µs per loop (mean ± std. dev. of 7 runs,1000 loops each)
In [14]: %%timeit
    ...: row_ind,col_ind = M.nonzero()
    ...: X = sparse.csr_matrix((M.shape[0]+2,M.shape[1]+2))
    ...: for i in [0,2]:
    ...:         if i!= 1 or j !=1:
    ...:             X += sparse.csr_matrix( (M.data,(
    ...: M.shape[0]+2,M.shape[1]+2))
    ...: X = X[1:-1,1:-1]
    ...: 
    ...: 
4.72 ms ± 92.4 µs per loop (mean ± std. dev. of 7 runs,100 loops each)
In [15]: timeit X = signal.convolve(M.A,mode='same')

85.9 µs ± 339 ns per loop (mean ± std. dev. of 7 runs,10000 loops each)