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

Question

格拉姆矩阵是结构X @ X.T对称的矩阵。在处理密集矩阵时，numpy.dot产品实现足够智能，可以识别自乘以利用对称性，从而加快计算速度（请参阅此）。scipy.sparse然而，使用矩阵时观察不到这样的效果：

\n\n

random.seed(0)\nX = random.randn(5,50)\nX[X < 1.5] = 0\nX = scipy.sparse.csr_matrix(X)\nprint(f\'sparsity of X: {100 * (1 - X.count_nonzero() / prod(X.shape)):5.2f} %\')\n# sparsity of X: 92.00 %\n\n%timeit X @ X.T\n# 248 \xc2\xb5s \xc2\xb1 10.8 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 1000 loops each)\n\nX2 = X.copy()\n%timeit X @ X2.T\n# 251 \xc2\xb5s \xc2\xb1 9.38 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 1000 loops each)\n

\n\n

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

\n\n

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

\n

Answer 1

感谢用户CJR的评论，我找到了一个令人满意的解决方案。事实上，我在 GitHub 上找到了一个库，它包装了 Python 的 MKL 例程mkl_sparse_spmm。该例程用于两个稀疏矩阵的快速乘法。所以我所要做的就是扩展库并为mkl_sparse_syrk. 而这正是我所做的。

\n\n

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

\n\n

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

\n\n

random.seed(0)\nX = random.randn(500, 5000)\nX[X < 0.8] = 0\nX = scipy.sparse.csr_matrix(X)\nprint(f\'X sparsity: {100 * (1 - X.count_nonzero() / prod(X.shape)):5.2f} %\')\n# X sparsity: 78.80 %\n\nexpected_result = (X @ X.T).toarray()\nexpected_result_triu = expected_result.copy()\nexpected_result_triu[tril_indices(expected_result.shape[0], k=-1)] = 0\n\nmkl_result1 = sparse_dot_mkl.dot_product_mkl(X, X.T)\nallclose(mkl_result1.toarray(), expected_result)\n# True\n\nmkl_result2 = sparse_dot_mkl.dot_product_transpose_mkl(X)\nallclose(mkl_result2.toarray(), expected_result_triu)\n# True\n\n%timeit X @ X.T\n# 197 ms \xc2\xb1 5.21 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n\n%timeit sparse_dot_mkl.dot_product_mkl(X, X.T)\n# 70.6 ms \xc2\xb1 593 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n\n%timeit sparse_dot_mkl.dot_product_transpose_mkl(X)\n# 34.2 ms \xc2\xb1 421 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n

\n\n

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

\n