python - 如何向量化生成 scipy 稀疏矩阵的代码？

Question

让我首先解释一下我想做什么。我正在尝试构建一个基于 m 包的推荐系统，每个包都有 n 个功能，存储在 m x n 稀疏矩阵 X 中。为此，我尝试运行 kNN 来获取包的 k 个最接近的匹配项。我想构建一个 m x m 稀疏矩阵 K，其中 K[i, j] 是行 X[i] 的点积如果X[j]是kNN为X[i]返回的包，则为和X[j]，否则为0。

这是我编写的代码:

X = ...
knn = NearestNeighbors(n_neighbors=self.n_neighbors, metric='l2')
knn.fit(X)
knn_indices = knn.kneighbors(X, return_distance=False)

m, k = X.shape[0], self.n_neighbors
K = lil_matrix((m, m))

for i, indices in enumerate(knn_indices):
    xi = X.getrow(i)
    for j in indices:
        xj = X.getrow(j)
        K[i, j] = xi.dot(xj.T)[0, 0]

我正在尝试找出如何提高效率。在我的场景中，m 约为 120 万，n 约为 50000，k 为 500，因此性能非常重要。

我填充K的最后一部分是我的程序的瓶颈。 getrow 看起来表现很差；根据 scipy 文档，它会复制行，因此 getrow 调用每次调用时最多可以复制 50k 个元素。另外，在最里面的循环中，我不知道如何获取 dot 的标量，而不是创建一个全新的 1x1 稀疏矩阵。

如何避免这些问题并加速/矢量化此代码的最后部分？谢谢。

Answer 1

In [21]: from scipy import sparse
In [22]: M = sparse.random(10,10,.2,'csr')
In [23]: M
Out[23]: 
<10x10 sparse matrix of type '<class 'numpy.float64'>'
    with 20 stored elements in Compressed Sparse Row format>

寻找M.A ，我选了这个小knn_indices用于测试的数组:

In [45]: knn = np.array([[4],[2],[],[1,3]])

你的双循环:

In [46]: for i, indices in enumerate(knn):
    ...:     xi = M[i,:]
    ...:     for j in indices:
    ...:         xj = M[j,:]
    ...:         print((xi*xj.T).A)
    ...:         
[[0.35494592]]
[[0.]]
[[0.08112133]]
[[0.56905781]]

内部循环可以压缩:

In [47]: for i, indices in enumerate(knn):
    ...:     xi = M[i,:]
    ...:     xj = M[indices,:]
    ...:     print((xi*xj.T).A)
    ...:         
[[0.35494592]]
[[0.]]
[]
[[0.08112133 0.56905781]]

并进行作业:

In [49]: k = sparse.lil_matrix((4,5))
In [50]: for i, indices in enumerate(knn):
    ...:     xi = M[i,:]
    ...:     for j in indices:
    ...:         xj = M[j,:]
    ...:         k[i,j] = (xi*xj.T)[0,0]
    ...:         
    ...:         
In [51]: k.A
Out[51]: 
array([[0.        , 0.        , 0.        , 0.        , 0.35494592],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.08112133, 0.        , 0.56905781, 0.        ]])

第二个循环

 k[i,indices] = (xi*xj.T)

做同样的事情。

也许可以用 i 做一些事情也循环，但这至少是一个开始。

那个knn不需要数组。由于内部列表长度不同，它无论如何都是一个对象数据类型。最好将其保留为列表。

填写此内容的替代方案 lil矩阵，将累加i , indices以及 coo 中的点积样式数组。

In [64]: r,c,d = [],[],[]
In [65]: for i, indices in enumerate(knn):
    ...:     xi = M[i,:]
    ...:     xj = M[indices,:]
    ...:     t = (xi*xj.T).data
    ...:     if len(t)>0:
    ...:         r.extend([i]*len(indices))
    ...:         c.extend(indices)
    ...:         d.extend(t)
    ...:         
In [66]: r,c,d
Out[66]: 
([0, 3, 3],
 [4, 1, 3],
 [0.3549459176547072, 0.08112132851228658, 0.5690578146292733])
In [67]: sparse.coo_matrix((d,(r,c))).A
Out[67]: 
array([[0.        , 0.        , 0.        , 0.        , 0.35494592],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.        , 0.08112133, 0.        , 0.56905781, 0.        ]])

在我的测试用例中，第二行没有任何非零值，需要在循环中进行额外的测试。我不知道这是否比lil快方法。