如何在密集的numpy矩阵和稀疏的scipy向量之间进行有效的矩阵乘法？

Question

使用@ 将密集的 numpy 矩阵与稀疏的 scipy 向量相乘是非常低效的。它似乎根本没有利用向量的稀疏性。

假设我们有

A = np.eye(5000)
x = np.eye(5000)[333]
x = scipy.sparse.coo_matrix(x).T # make it a sparse vector

然后使用@ 进行乘法运算：

%timeit A @ x
8 ms ± 78.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

让我们自己写一个非常糟糕的稀疏乘法：

def mult_dense_with_sparse(A, x):
  return (A[:,x.nonzero()[0]] @ x.toarray()[x.nonzero()[0]]).T[0]

你瞧：

%timeit mult_dense_with_sparse(A, x)
50.3 µs ± 45.3 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

它快得多！即使这个实现首先通过再次添加所有零来创建一个密集向量，然后再次删除所有零......所以我想知道，如果不是@，我怎样才能将密集的 numpy 矩阵与稀疏的 scipy 相乘矢量有效？这样的基本操作肯定是 scipy 的一部分吗？

编辑：在其他问题中提供的解决方案没有帮助，因为它与@ 一样低效：

%timeit scipy.sparse.csr_matrix.dot(A, x)
7.97 ms ± 113 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

在响应 Hameer Abbasi 中编辑 2：

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   101                                                       @profile
   102                                                       def ratio(self, n):
   103        80         51.0      0.6      0.0                  s = n + 1
   104        80   11401854.0 142523.2     16.1                  self.sfpc = self.scolPCA(s) # sparse first principal component
   106        80     351898.0   4398.7      0.5                  wSums = (self.signals[:,self.sfpc.nonzero()[0]] @ self.sfpc.toarray()[self.sfpc.nonzero()[0]]).T[0]
   108        80   56487433.0 706092.9     79.7                  wSums = self.sfpc.T.dot(self.signals.T)[0]
   110        80    2521189.0  31514.9      3.6                  self.Flag, self.threshold, self.incline, self.deltaVar = self.actFunctOpt(list(wSums))
   111        80        160.0      2.0      0.0                  return  self.deltaVar / (2 + s)

在这里您可以看到，我们的“hack”大约占用了该函数 0.5% 的时间，而使用 dot 占用了该函数 79.7% 的时间。

Answer 1

在您的示例中，A 的类型为 np.ndarray，x 的类型为 scipy.sparse.coo_matrix。

如果您正在寻找最简单的答案来加快速度，这里是：

import numpy as np
import scipy.sparse as sps
A = np.random.rand(5000, 5000)
v = np.zeros((5000, 1))
v[333] = 1
v_s = sps.csc_matrix(v)
%timeit v_s.T.dot(A).T # 46.6 µs ± 1.11 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit A @ v # 6.75 ms ± 29.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

但是，如果您想了解此答案背后的机制：目前，由于 bug，运算符 @ 不支持 NumPy 数组的重载。此外，scipy.sparse.coo_matrix 甚至没有尝试覆盖@，scipy.sparse.csr_matrix 的矩阵乘法更快（无论如何，coo_matrix 将首先转换为csr_matrix 进行乘法运算）。

NumPy 将您的稀疏向量转换为 NumPy 数组，然后执行密集乘法。

但是，csr_matrix.dot 存在并支持与密集 NumPy 数组的乘法。所以我们使用属性A @ B = (B.T @ A.T).T，以及csc_matrix.T产生csr_matrix的事实来产生上述优化代码。