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

格拉姆矩阵是结构矩阵X @ X.T这当然是对称的。当处理稠密矩阵时，numpy.dot产品实现足够智能，可以识别自乘以利用对称性，从而加快计算速度（请参阅this https://stackoverflow.com/a/50734430/1444073）。然而，使用时观察不到这样的效果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中计算稀疏格拉姆矩阵的最快方法是什么？值得注意的是，仅计算下三角形（或等效的上三角形）就足够了。

我读过多次，使用天际线格式 https://en.wikipedia.org/wiki/Skyline_matrix对于对称矩阵非常有效，但是 scipy 不支持天际线格式。相反，人们把矛头指向pysparse http://pysparse.sourceforge.net/很多次，但似乎pysparse已经停产很久了，并且不支持Python 3。至少，我的Anaconda由于与Python 3的兼容性问题而拒绝安装pysparse。

感谢用户CJR的评论，我找到了一个令人满意的解决方案。事实上，我发现GitHub 上的一个库 https://github.com/flatironinstitute/sparse_dot它包装了 MKL 例程mkl_sparse_spmm对于Python。该例程用于两个稀疏矩阵的快速乘法。所以我所要做的就是扩展库并提供类似的包装器mkl_sparse_syrk。这正是我做了什么 https://github.com/kostrykin/sparse_dot.

我还需要添加一些评论，之后我将向原始项目提交拉取请求。

然而，以下是性能结果，令人印象深刻：

random.seed(0)
X = random.randn(500, 5000)
X[X < 0.8] = 0
X = scipy.sparse.csr_matrix(X)
print(f'X sparsity: {100 * (1 - X.count_nonzero() / prod(X.shape)):5.2f} %')
# X sparsity: 78.80 %

expected_result = (X @ X.T).toarray()
expected_result_triu = expected_result.copy()
expected_result_triu[tril_indices(expected_result.shape[0], k=-1)] = 0

mkl_result1 = sparse_dot_mkl.dot_product_mkl(X, X.T)
allclose(mkl_result1.toarray(), expected_result)
# True

mkl_result2 = sparse_dot_mkl.dot_product_transpose_mkl(X)
allclose(mkl_result2.toarray(), expected_result_triu)
# True

%timeit X @ X.T
# 197 ms ± 5.21 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit sparse_dot_mkl.dot_product_mkl(X, X.T)
# 70.6 ms ± 593 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit sparse_dot_mkl.dot_product_transpose_mkl(X)
# 34.2 ms ± 421 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

使用 MKL 中的通用点积代替 scipy 中的点积实现会产生加速 279%。使用专门的产品进行 Gram 矩阵计算可以得到加速 576%。这是巨大的。