Apply a convolution to a scipy.sparse matrix

Question

I try to compute a convolution on a scipy.sparse matrix. Here is the code:

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

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

Which produces the following error:

ValueError: volume and kernel should have the same dimensionality

Computing scipy.signal.convolve(M.todense(), kernel, mode='same') provides the expected result. However, I would like to keep the computation sparse.

More generally speaking, my goal is to compute the 1-hop neighbourhood sum of the sparse matrix M. If you have any good idea how to calculate this on a sparse matrix, I would love to hear it !

EDIT:

I just tried a solution for this specific kernel (sum of neighbors) that is not really faster than the dense version (I didn't try in a very high dimension though). Here is the code:

row_ind, col_ind = M.nonzero() 
X = scipy.sparse.csr_matrix((M.shape[0]+2, M.shape[1]+2))
for i in [0, 1, 2]:
    for j in [0, 1, 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]

Answer 1

In [1]: from scipy import sparse, signal
In [2]: M = sparse.csr_matrix([[0,1,0,0],[1,0,0,1],[1,0,1,0],[0,0,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.],
       [1., 3., 1., 2.],
       [1., 2., 1., 1.]])

Why do posters show runnable code, but not the results? Most of us can't run code like this in our heads.

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

Your alternative - while the result is a sparse matrix, all values are filled. Even if M is larger and sparser, X will be denser.

In [7]: row_ind, col_ind = M.nonzero()
   ...: X = sparse.csr_matrix((M.shape[0]+2, M.shape[1]+2))
   ...: for i in [0, 1, 2]:
   ...:     for j in [0, 1, 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., 2., 1.],
       [2., 4., 3., 1.],
       [1., 3., 1., 2.],
       [1., 2., 1., 1.]])

Here's an alternative that builds the coo style inputs, and only makes the matrix at the end. Keep in mind that repeated coordinates are summed. That's handy in FEM stiffness matrix construction, and fits nicely here as well.

In [10]: row_ind, col_ind = M.nonzero()
    ...: data, row, col = [],[],[]
    ...: for i in [0, 1, 2]:
    ...:     for j in [0, 1, 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, 1, 2, 1],
       [2, 4, 3, 1],
       [1, 3, 1, 2],
       [1, 2, 1, 1]])

===

My approach is noticeably faster (but still well behind the dense convolution). sparse.csr_matrix(...) is pretty slow, so it isn't a good idea to do repeatedly. And sparse addition isn't very good either.

In [13]: %%timeit
    ...: row_ind, col_ind = M.nonzero()
    ...: data, row, col = [],[],[]
    ...: for i in [0, 1, 2]:
    ...:     for j in [0, 1, 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]
    ...: 
    ...: 
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, 1, 2]:
    ...:     for j in [0, 1, 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]
    ...: 
    ...: 
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, kernel, mode='same')

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