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

问题描述

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，但似乎很久以前就已经停产了，并且没有对Python3的支持。至少，我的Anaconda因为与Python3的兼容性问题而拒绝安装pyparse。

---

解决方案

感谢cjr用户的评论，我得出了一个令人满意的解决方案。事实上，我找到了a library on GitHub，它包装了用于Python的MKL例程mkl_sparse_spmm。此例程用于两个稀疏矩阵的快速相乘。所以我所要做的就是扩展这个库，并为mkl_sparse_syrk提供一个类似的包装器。而这正是what I did。

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

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

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%。这是件大事。