稀疏矩阵上的逐行外积答案

【问题标题】：Row-wise outer product on sparse matrices稀疏矩阵上的逐行外积
【发布时间】：2019-11-27 16:43:10
【问题描述】：

给定两个稀疏的 scipy 矩阵 A, B 我想计算逐行外积。

我可以通过多种方式使用 numpy 做到这一点。最简单的可能是

np.einsum('ij,ik->ijk', A, B).reshape(n, -1)

或

(A[:, :, np.newaxis] * B[:, np.newaxis, :]).reshape(n, -1)

其中n 是A 和B 中的行数。

然而，就我而言，通过密集矩阵会消耗太多 RAM。因此，我发现的唯一选择是使用 python 循环：

sp.sparse.vstack((ra.T@rb).reshape(1,-1) for ra, rb in zip(A,B)).tocsr()

虽然使用较少的 RAM，但速度非常慢。

因此，我的问题是，是否有一种稀疏（RAM 高效）方法来获取两个矩阵的逐行外积，从而使事物保持矢量化？

（类似的问题是numpy elementwise outer product with sparse matrices，但那里的所有答案都经过密集矩阵。）

【问题讨论】：

A 和 B 的典型尺寸是多少？
@Divakar 我有 1000,000 行和 768 列。每列的典型非零值约为 25。
遍历稀疏矩阵的行很慢，因为它需要创建一个新的 (1,n) 稀疏矩阵。 @ 应该足够快，重塑应该没问题，因为它简化了 csr indptr。 vstack 结合了所有输入矩阵的 coo 属性，对于这么多行可能会很慢。
lil 格式有一个快速行 view 方法，但它仍然必须转换为 csr 格式才能进行乘法运算。
所以你期望一个有 1000,000 行和 768*768 列的矩阵？

标签： python numpy matrix scipy sparse-matrix

【解决方案1】：

我们可以直接计算结果的csr表示。它不是超快的（在 100,000x768 上约 3 秒）但可能没问题，具体取决于您的用例：

import numpy as np
import itertools
from scipy import sparse

def spouter(A,B):
    N,L = A.shape
    N,K = B.shape
    drows = zip(*(np.split(x.data,x.indptr[1:-1]) for x in (A,B)))
    data = [np.outer(a,b).ravel() for a,b in drows]
    irows = zip(*(np.split(x.indices,x.indptr[1:-1]) for x in (A,B)))
    indices = [np.ravel_multi_index(np.ix_(a,b),(L,K)).ravel() for a,b in irows]
    indptr = np.fromiter(itertools.chain((0,),map(len,indices)),int).cumsum()
    return sparse.csr_matrix((np.concatenate(data),np.concatenate(indices),indptr),(N,L*K))

A = sparse.random(100,768,0.03).tocsr()
B = sparse.random(100,768,0.03).tocsr()

print(np.all(np.einsum('ij,ik->ijk',A.A,B.A).reshape(100,-1) == spouter(A,B).A))

A = sparse.random(100000,768,0.03).tocsr()
B = sparse.random(100000,768,0.03).tocsr()

from time import time
T = time()
C = spouter(A,B)
print(time()-T)

示例运行：

True
3.1073222160339355

【讨论】：

这似乎比我的版本快。你介意解释一下它是如何工作的吗？
它逐行工作。在 A[i] 外部 B[i] 中找到非零值意味着将 A[i] 和 B[i] 中的所有非零值配对 -> value = A[i].data[j]*B[i].data[k]在 (A[i].indices[j], B[i].indices[k]);然后必须将索引对转换为线性索引。
我不明白为什么它比我的版本 (3) 快得多，因为它们都对所有行进行了完整的 python 循环？
只要最里面的循环是矢量化的，循环不是很好但也不可怕。因此，差异归结为最内层循环的实现方式。您正在使用 matmul，这是一个非常通用的函数，而我使用的代码只能执行这种特殊情况（两个（稀疏）向量的外积）。另外我认为你的循环为每一行创建了完整的稀疏矩阵对象，我认为这是相当多的开销