向量和 csr_matrix 的左乘

Question

我有一个大小为 50*10000 的 csr 格式稀疏矩阵 A 和一个大小为 50 的 np.array v，我想为此计算乘积 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)
```

1. 转置 v 并使用乘法和求和法

我正在使用csr_matrix.multiply() 方法，它执行逐点乘法。我会总结这些行。

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

1. 将 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 返回一个新矩阵，但采用 csc 格式。否则变化很小。但它不是像密集转置那样的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)

Sparse 针对非常稀疏的大型矩阵的矩阵乘法进行了优化。大多数其他操作使用密集数组会更快（只要它们适合内存）。

稀疏的默认乘法是矩阵

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__。稀疏的dot 也可以。

* 和dot 时间相同，@ 慢一点。 在 [758] 中：timeit A.T*(v)
每个循环 95.6 µs ± 424 ns（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

@ 运算符也适用于 scipy 稀疏矩阵，与您的方法的简短和粗略的性能比较表明 @ 运算符的执行类似于您的方法 1