Объединение разреженных матриц Scipy без медленного scipy.sparse.vstack

Question

Мне нужно вычислить массив размером 1xN в пределах al oop и сложить каждый новый массив поверх предыдущего. Длина l oop равна 1, M, и функция работает аналогично следующему:

import numpy as np
import scipy.sparse as sp

mat = np.random.uniform(size=(1,N))
mat_sp = sp.coo_matrix(mat)

for i in range(1,M):

    mat_new = np.random.uniform(size=(1,N))
    mat_sp_new = sp.coo_matrix(mat_new)

    mat_sp = sp.vstack((mat_sp,mat_sp_new))

В результате получается матрица MxN. Однако выполнение этой операции с scipy.sparse.vstack выполняется очень медленно. Сравнение с numpy .vstack или даже просто с предварительно выделенной матрицей и N = 10000:

In [1]: %timeit sp.vstack((mat_sp,mat_sp)) 
315 µs ± 7.89 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [2]: %timeit np.vstack((mat,mat))
8.63 µs ± 87.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [3]: mat_final = np.zeros((2,10000))   

In [4]: %timeit mat_final[0]=mat; mat_final[1]=mat                               
4.03 µs ± 92.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Моя проблема в том, что иногда мне понадобится моя конечная матрица, чтобы иметь размеры до M = 10 ⁵ и N = 10 ⁶ , что с NumPy даст MemoryError.

Есть ли способ наложения векторов разреженных строк, который не не так ли медленно?

Answer 1

Использование повторяющегося stack, как в mat_sp = sp.vstack((mat_sp,mat_sp_new)), является плохой идеей, даже если используется np.vstack. Гораздо лучше собрать все массивы в списке и использовать stack только один раз. Добавление списка намного эффективнее, когда используется итеративно, а stacks предназначен для работы со всем списком объектов, а не только с двумя одновременно. Список приложений работает на месте со ссылками на объекты; операции с массивами каждый раз возвращают новый массив.

Позвольте мне проиллюстрировать:

In [1]: from scipy import sparse 

функция для создания разреженного 1d массива:

In [2]: def foo(): 
   ...:     x = np.zeros(10,int) 
   ...:     idx = np.random.randint(0,10,size=4) 
   ...:     x[idx] = idx 
   ...:     return x 
   ...:                                                                                        
In [3]: foo()                                                                                  
Out[3]: array([0, 1, 2, 3, 0, 0, 0, 0, 0, 9])
In [4]: foo()                                                                                  
Out[4]: array([0, 0, 0, 3, 0, 0, 0, 0, 0, 9])

Добавление простого списка :

In [5]: alist = []                                                                             
In [6]: for i in range(4): 
   ...:     alist.append(foo()) 
   ...:                                                                                        
In [7]: alist                                                                                  
Out[7]: 
[array([0, 0, 2, 0, 0, 0, 0, 0, 8, 9]),
 array([0, 0, 0, 3, 0, 5, 0, 0, 0, 9]),
 array([0, 1, 0, 0, 0, 0, 6, 7, 0, 0]),
 array([0, 1, 0, 0, 0, 0, 0, 0, 8, 0])]

плотный массив из этого:

In [8]: np.vstack(alist)                                                                       
Out[8]: 
array([[0, 0, 2, 0, 0, 0, 0, 0, 8, 9],
       [0, 0, 0, 3, 0, 5, 0, 0, 0, 9],
       [0, 1, 0, 0, 0, 0, 6, 7, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0, 8, 0]])

разреженная матрица из этого массива:

In [9]: M = sparse.coo_matrix(_)                                                               
In [10]: M                                                                                     
Out[10]: 
<4x10 sparse matrix of type '<class 'numpy.int64'>'
    with 11 stored elements in COOrdinate format>
In [11]: M                                                                                     
Out[11]: 
<4x10 sparse matrix of type '<class 'numpy.int64'>'
    with 11 stored elements in COOrdinate format>
In [12]: print(M)                                                                              
  (0, 2)    2
  (0, 8)    8
  (0, 9)    9
  (1, 3)    3
  (1, 5)    5
  (1, 9)    9
  (2, 1)    1
  (2, 6)    6
  (2, 7)    7
  (3, 1)    1
  (3, 8)    8

Альтернатива - создать разреженную матрицу из каждого массива и объединить те:

In [13]: alist = []                                                                            
In [14]: for i in range(4): 
    ...:     alist.append(sparse.coo_matrix(foo())) 
    ...:                                                                                       
In [15]: alist                                                                                 
Out[15]: 
[<1x10 sparse matrix of type '<class 'numpy.int64'>'
    with 2 stored elements in COOrdinate format>,
 <1x10 sparse matrix of type '<class 'numpy.int64'>'
    with 3 stored elements in COOrdinate format>,
 <1x10 sparse matrix of type '<class 'numpy.int64'>'
    with 2 stored elements in COOrdinate format>,
 <1x10 sparse matrix of type '<class 'numpy.int64'>'
    with 3 stored elements in COOrdinate format>]
In [16]: M1=sparse.vstack(alist)                                                               
In [17]: M1                                                                                    
Out[17]: 
<4x10 sparse matrix of type '<class 'numpy.longlong'>'
    with 10 stored elements in COOrdinate format>

Разреженная матрица (формат coo) хранит свою информацию в 3 массивах (ср. с [12]):

In [18]: M.data                                                                                
Out[18]: array([2, 8, 9, 3, 5, 9, 1, 6, 7, 1, 8])
In [19]: M.row                                                                                 
Out[19]: array([0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3], dtype=int32)
In [20]: M.col                                                                                 
Out[20]: array([2, 8, 9, 3, 5, 9, 1, 6, 7, 1, 8], dtype=int32)

Обратите внимание, как первые 3 значения в эти массивы соответствуют ненулевым элементам:

array([0, 0, 2, 0, 0, 0, 0, 0, 8, 9]

In [21]: m1 =sparse.coo_matrix(np.array([0, 0, 2, 0, 0, 0, 0, 0, 8, 9])) 
In [24]: m1.data                                                                               
Out[24]: array([2, 8, 9])
In [25]: m1.row                                                                                
Out[25]: array([0, 0, 0], dtype=int32)
In [26]: m1.col                                                                                
Out[26]: array([2, 8, 9], dtype=int32)

Распространенным способом построения разреженной матрицы является прямой сбор массивов data, row, col, возможно, в виде списков, а затем подача к coo_matrix:

sparse.coo_matrix((data, (row, col)), shape=(N,M))

sparse.vstack передает задачу на sparse.bmat. Посмотрите на код bmat, чтобы увидеть, как он собирает атрибуты data/row/col для каждого из входных данных в составные массивы, а затем создает матрицу coo.

Генератор образцов матрицы:

In [38]: data, row, col = [], [], []                                                           
In [39]: for i in range(4): 
    ...:     idx = np.random.randint(0,10,size=4) 
    ...:     data.append(idx) 
    ...:     row.append(np.ones(len(idx),int)*i) 
    ...:     col.append(idx) 
    ...: data = np.hstack(data) 
    ...: row = np.hstack(row) 
    ...: col = np.hstack(col) 
    ...: M = sparse.coo_matrix((data, (row, col)), shape=(4,10)) 
    ...:  
    ...:                                                                                       
In [40]: M                                                                                     
Out[40]: 
<4x10 sparse matrix of type '<class 'numpy.int64'>'
    with 16 stored elements in COOrdinate format>
In [41]: data                                                                                  
Out[41]: array([5, 8, 6, 2, 6, 3, 5, 8, 6, 1, 3, 5, 0, 2, 7, 0])
In [42]: row                                                                                   
Out[42]: array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3])
In [43]: M.A                                                                                   
Out[43]: 
array([[0, 0, 2, 0, 0, 5, 6, 0, 8, 0],
       [0, 0, 0, 3, 0, 5, 6, 0, 8, 0],
       [0, 1, 0, 3, 0, 5, 6, 0, 0, 0],
       [0, 0, 2, 0, 0, 0, 0, 7, 0, 0]])