在Python中计算稀疏Gram矩阵最快的方法是什么？

Question

Gram矩阵是结构X @ X.T的矩阵，它当然是对称的。在处理密集矩阵时，numpy.dot乘积实现足够智能，可以识别自乘以利用对称性，从而加快计算速度（请参见this）。但是，当使用scipy.sparse矩阵时，无法观察到这种效果：

random.seed(0)
X = random.randn(5,50)
X[X < 1.5] = 0
X = scipy.sparse.csr_matrix(X)
print(f'sparsity of X: {100 * (1 - X.count_nonzero() / prod(X.shape)):5.2f} %')
# sparsity of X: 92.00 %

%timeit X @ X.T
# 248 µs ± 10.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

X2 = X.copy()
%timeit X @ X2.T
# 251 µs ± 9.38 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

[因此，我想知道：在Python中计算稀疏Gram矩阵的最快方法是什么？值得注意的是，仅计算下部（或等效的上部）三角形就足够了。

我已经读过很多遍了，使用skyline format对于对称矩阵非常有效，但是scipy不支持天际线格式。取而代之的是，人们多次指向pysparse，但是pysparse似乎早就停产了，并且不支持Python3。至少，由于与Python 3的兼容性问题，我的Anaconda拒绝安装pysparse。

Answer 1

通过用户CJR的评论，我制定了一个令人满意的解决方案。实际上，我发现a library on GitHub包装了Python的MKL例程mkl_sparse_spmm。此例程用于两个矩阵的快速乘法。因此，我要做的就是扩展库并为mkl_sparse_syrk提供类似的包装器。这就是what I did。