Блог

Я пытаюсь вычислить свертку на матрице scipy.sparse. Вот код:

 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')
  

Что приводит к следующей ошибке:

 ValueError: volume and kernel should have the same dimensionality
  

Вычисления scipy.signal.convolve(M.todense(), kernel, mode='same') дают ожидаемый результат. Тем не менее, я хотел бы сохранить вычисление разреженным.

В более общем плане, моя цель — вычислить сумму окрестностей за 1 переход разреженной матрицы M. Если у вас есть хорошая идея, как вычислить это на разреженной матрице, я хотел бы это услышать!

Редактировать:

Я только что попробовал решение для этого конкретного ядра (сумма соседей), которое на самом деле не быстрее, чем плотная версия (хотя я не пробовал в очень высоком измерении). Вот код:

 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]
  

 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.]])
  

Почему на плакатах отображается выполняемый код, но не результаты? Большинство из нас не могут запускать подобный код в своих головах.

 In [5]: M.A
Out[5]: 
array([[0, 1, 0, 0],
       [1, 0, 0, 1],
       [1, 0, 1, 0],
       [0, 0, 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, 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.]])
  

Вот альтернатива, которая создает входные coo данные стиля и создает матрицу только в конце. Имейте в виду, что повторяющиеся координаты суммируются. Это удобно при построении матрицы жесткости FEM и хорошо подходит и здесь.

 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]])
  

===

Мой подход заметно быстрее (но все еще значительно отстает от плотной свертки). sparse.csr_matrix(...) работает довольно медленно, поэтому повторять это не рекомендуется. И разреженное добавление тоже не очень хорошо.

 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)