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]

`sparse` matrices are not a subclass of `ndarray`, so `numpy` and other `scipy` modules don't generally handle them correctly. That's why you need to use `todense`. — hpaulj, Oct 21 '20 at 16:28
Sparse matrices are best for matrix multiplication and math that does not change sparsity. Adding matrices is relatively slow. So is repeatedly creating a matrix, as you do in the loop. But you may instead be able to collect all those `M.data`, `row_ind+i` etc values in `coo` style arrays, and do one matrix construction at the end. Duplicate values are summed. — hpaulj, Oct 21 '20 at 16:34

score 1 · Accepted Answer · answered Oct 21 '20 at 17:09

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)

Right, calling the sparse.csr_matrix constructor only once is far better than my trivial solution ! I guess it is the best solution given this specific kernel. If M is big (and sparse), then this solution is way faster than the dense version (using convolve) as well. — Robin, Oct 22 '20 at 09:42

Apply a convolution to a scipy.sparse matrix

1 Answers1