Левое Умножение вектора и csr_matrix

Question

У меня есть разреженная матрица формата csr A размера 50 * 10000 и массив np.array v размера 50, для которого я хотел бы рассчитать произведение v.dot(A). Как мне сделать это эффективно?

Конечно, запуск v.dot(A) не очень хорошая идея, потому что scipy.sparse не поддерживает операции np. К сожалению, насколько мне известно, scipy.sparse не имеет функции для умножения влево матрицы на вектор.

Я пробовал следующие методы, но все они кажутся довольно дорогостоящими:

1. Транспонировать A и использовать стандарт .dot

Я переставляю A и затем использую метод .dot. Это умножает A.T на v в качестве вектора столбца.

```
>>> A = sparse.csr_matrix([[1, 2, 0, 0],
...                        [0, 0, 3, 4]])
>>> v = np.array([1, 1])
>>> A.T.dot(v)
array([1, 2, 3, 4], dtype=int32)
```

Транспонирование v и использование методов умножения и суммирования

Я использую метод csr_matrix.multiply(), который выполняет точечное умножение. Я буду суммировать по строкам.

>>> vt = v[np.newaxis].T
>>> A.multiply(vt).sum(axis=0)
matrix([[1, 2, 3, 4]], dtype=int32)

Превращение v в разреженную матрицу и использование метода .dot

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

>>> sparse_v = sparse.csr_matrix(v)
>>> sparse_v.dot(A).todense()
matrix([[1, 2, 3, 4]], dtype=int32)

Метод 1, безусловно, самый быстрый, но метод .T все еще очень трудоемкий. Нет ли лучшего способа выполнить левое умножение на разреженные матрицы?

Answer 1

In [746]: A = sparse.csr_matrix([[1, 2, 0, 0], 
     ...: ...                        [0, 0, 3, 4]])                                                    
In [747]: A                                                                                            
Out[747]: 
<2x4 sparse matrix of type '<class 'numpy.longlong'>'
    with 4 stored elements in Compressed Sparse Row format>
In [748]: print(A)                                                                                     
  (0, 0)    1
  (0, 1)    2
  (1, 2)    3
  (1, 3)    4
In [749]: v = np.array([1, 1])                                                                         

A.T возвращает новую матрицу, но в формате cs c. В противном случае изменение минимально. Но это не view как плотная транспонирование.

In [750]: A.T                                                                                          
Out[750]: 
<4x2 sparse matrix of type '<class 'numpy.longlong'>'
    with 4 stored elements in Compressed Sparse Column format>
In [755]: timeit A.T                                                                                   
82.9 µs ± 542 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

На самом деле время довольно хорошее по сравнению со временем, необходимым для создания исходной разреженной матрицы:

In [759]: timeit A = sparse.csr_matrix([[1, 2, 0, 0],   [0, 0, 3, 4]])                                 
349 µs ± 356 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Разреженная оптимизирована для умножения матриц больших очень разреженных матриц. Большинство других операций выполняются быстрее с плотными массивами (при условии, что они помещаются в памяти).

Умножением по умолчанию для разреженных является матрица

In [752]: A.T*v                                                                                        
Out[752]: array([1, 2, 3, 4], dtype=int64)
In [753]: A.T.dot(v)                                                                                   
Out[753]: array([1, 2, 3, 4], dtype=int64)
In [754]: A.T@(v)                                                                                      
Out[754]: array([1, 2, 3, 4], dtype=int64)

@, __matmul__ делегаты __mul__ , Sparse dot делает то же самое.

* и dot время одинаковые, @ немного медленнее. В [758]: время AT * (v)
95,6 мкс ± 424 нс на л oop (среднее ± стандартное отклонение из 7 циклов, 10000 циклов каждый)

, удаляя транспонирование из время:

In [760]: %%timeit A1=A.T 
     ...: A1*v                                                                                              
7.77 µs ± 22.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Ваш второй подход использует разреженное поэлементное умножение:

In [774]: A.multiply(v[:,None]).sum(axis=0)                                                            
Out[774]: matrix([[1, 2, 3, 4]], dtype=int64)

Это не так эффективно, как умножение матриц

In [775]: A.multiply(v[:,None])                                                                        
Out[775]: 
<2x4 sparse matrix of type '<class 'numpy.longlong'>'
    with 4 stored elements in COOrdinate format>
In [776]: timeit A.multiply(v[:,None])                                                                 
147 µs ± 1.14 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

В результате получается еще одна разреженная матрица , sum(axis=0) фактически реализовано как матричное умножение. Матрица extractor для этой суммы равна sparse.csr_matrix([1,1]). Но это всего лишь sparse_v в вашем последнем примере.

In [787]: A.sum(axis=0)                                                                                
Out[787]: matrix([[1, 2, 3, 4]], dtype=int64)
In [788]: timeit A.sum(axis=0)                                                                         
242 µs ± 191 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Все эти моменты времени следует рассматривать с осторожностью. A маленький и не очень редкий. Сравните эту сумму строки для гораздо большей разреженной матрицы:

In [799]: Ab = sparse.random(1000,1000,.001,'csr')                                                     
In [800]: timeit Ab.sum(axis=0)                                                                        
269 µs ± 12.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [801]: timeit Ab.T*np.ones(1000)                                                                    
118 µs ± 216 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Answer 2

Оператор @ также работает для скудных разреженных матриц, а краткое и грубое сравнение производительности с вашими методами показывает, что оператор @ работает аналогично вашему методу 1