问题描述
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)